diff --git a/python/core/raster/qgsrasterblock.sip b/python/core/raster/qgsrasterblock.sip index 15546f5d7c00..75a620adaa6f 100644 --- a/python/core/raster/qgsrasterblock.sip +++ b/python/core/raster/qgsrasterblock.sip @@ -229,6 +229,10 @@ class QgsRasterBlock void applyNoDataValues( const QgsRasterRangeList & rangeList ); + /** apply band scale and offset to raster block values + * @@note added in 2.3 */ + void applyScaleOffset( double scale, double offset ); + /** \brief Get error */ QgsError error() const; diff --git a/python/core/raster/qgsrasterdataprovider.sip b/python/core/raster/qgsrasterdataprovider.sip index 55ceb5e5a9dc..b4c1f4b3b853 100644 --- a/python/core/raster/qgsrasterdataprovider.sip +++ b/python/core/raster/qgsrasterdataprovider.sip @@ -56,6 +56,14 @@ class QgsRasterDataProvider : QgsDataProvider, QgsRasterInterface virtual QString colorInterpretationName( int theBandNo ) const; + /** Read band scale for raster value + * @@note added in 2.3 */ + virtual double bandScale( int bandNo ) const; + + /** Read band offset for raster value + * @@note added in 2.3 */ + virtual double bandOffset( int bandNo ) const; + /** Get block size */ virtual int xBlockSize() const; virtual int yBlockSize() const; diff --git a/src/core/raster/qgsrasterblock.cpp b/src/core/raster/qgsrasterblock.cpp index 08127101f23c..a9eeffd15d06 100644 --- a/src/core/raster/qgsrasterblock.cpp +++ b/src/core/raster/qgsrasterblock.cpp @@ -706,6 +706,20 @@ bool QgsRasterBlock::convert( QGis::DataType destDataType ) return true; } +void QgsRasterBlock::applyScaleOffset( double scale, double offset ) +{ + if ( isEmpty() ) return; + if ( !typeIsNumeric( mDataType ) ) return; + if ( scale == 1.0 && offset == 0.0 ) return; + + qgssize size = (qgssize) mWidth * mHeight; + for ( qgssize i = 0; i < size; ++i ) + { + if ( !isNoData( i ) ) setValue( i, value( i ) * scale + offset ); + } + return; +} + void QgsRasterBlock::applyNoDataValues( const QgsRasterRangeList & rangeList ) { if ( rangeList.isEmpty() ) diff --git a/src/core/raster/qgsrasterblock.h b/src/core/raster/qgsrasterblock.h index b046a25c3304..d720d85cb4d7 100644 --- a/src/core/raster/qgsrasterblock.h +++ b/src/core/raster/qgsrasterblock.h @@ -291,6 +291,10 @@ class CORE_EXPORT QgsRasterBlock void applyNoDataValues( const QgsRasterRangeList & rangeList ); + /** apply band scale and offset to raster block values + * @@note added in 2.3 */ + void applyScaleOffset( double scale, double offset ); + /** \brief Get error */ QgsError error() const { return mError; } diff --git a/src/core/raster/qgsrasterdataprovider.cpp b/src/core/raster/qgsrasterdataprovider.cpp index 345d4f31f276..70a404460ac9 100644 --- a/src/core/raster/qgsrasterdataprovider.cpp +++ b/src/core/raster/qgsrasterdataprovider.cpp @@ -208,6 +208,8 @@ QgsRasterBlock * QgsRasterDataProvider::block( int theBandNo, QgsRectangle cons readBlock( theBandNo, theExtent, theWidth, theHeight, block->bits() ); } + // apply scale and offset + block->applyScaleOffset( bandScale( theBandNo ), bandOffset( theBandNo ) ); // apply user no data values block->applyNoDataValues( userNoDataValues( theBandNo ) ); return block; diff --git a/src/core/raster/qgsrasterdataprovider.h b/src/core/raster/qgsrasterdataprovider.h index 21479bea04a1..e0f6f5c37e32 100644 --- a/src/core/raster/qgsrasterdataprovider.h +++ b/src/core/raster/qgsrasterdataprovider.h @@ -161,6 +161,13 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider, public QgsRast return colorName( colorInterpretation( theBandNo ) ); } + /** Read band scale for raster value + * @@note added in 2.3 */ + virtual double bandScale( int bandNo ) const { Q_UNUSED( bandNo ); return 1.0; } + /** Read band offset for raster value + * @@note added in 2.3 */ + virtual double bandOffset( int bandNo ) const { Q_UNUSED( bandNo ); return 0.0; } + // TODO: remove or make protected all readBlock working with void* /** Read block of data using given extent and size. */ diff --git a/src/providers/gdal/qgsgdalprovider.cpp b/src/providers/gdal/qgsgdalprovider.cpp index c40144272771..5e2dd7c9ea3d 100644 --- a/src/providers/gdal/qgsgdalprovider.cpp +++ b/src/providers/gdal/qgsgdalprovider.cpp @@ -387,6 +387,8 @@ QgsRasterBlock* QgsGdalProvider::block( int theBandNo, const QgsRectangle &theEx block->setIsNoDataExcept( subRect ); } readBlock( theBandNo, theExtent, theWidth, theHeight, block->bits() ); + // apply scale and offset + block->applyScaleOffset( bandScale( theBandNo ), bandOffset( theBandNo ) ); block->applyNoDataValues( userNoDataValues( theBandNo ) ); return block; } @@ -1067,55 +1069,45 @@ int QgsGdalProvider::capabilities() const return capability; } -QGis::DataType QgsGdalProvider::dataTypeFormGdal( int theGdalDataType ) const -{ - switch ( theGdalDataType ) - { - case GDT_Unknown: - return QGis::UnknownDataType; - break; - case GDT_Byte: - return QGis::Byte; - break; - case GDT_UInt16: - return QGis::UInt16; - break; - case GDT_Int16: - return QGis::Int16; - break; - case GDT_UInt32: - return QGis::UInt32; - break; - case GDT_Int32: - return QGis::Int32; - break; - case GDT_Float32: - return QGis::Float32; - break; - case GDT_Float64: - return QGis::Float64; - break; - case GDT_CInt16: - return QGis::CInt16; - break; - case GDT_CInt32: - return QGis::CInt32; - break; - case GDT_CFloat32: - return QGis::CFloat32; - break; - case GDT_CFloat64: - return QGis::CFloat64; - break; - } - return QGis::UnknownDataType; -} - QGis::DataType QgsGdalProvider::srcDataType( int bandNo ) const { GDALRasterBandH myGdalBand = GDALGetRasterBand( mGdalDataset, bandNo ); GDALDataType myGdalDataType = GDALGetRasterDataType( myGdalBand ); - return dataTypeFromGdal( myGdalDataType ); + QGis::DataType myDataType = dataTypeFromGdal( myGdalDataType ); + + // define if the band has scale and offset to apply + double myScale = bandScale( bandNo ); + double myOffset = bandOffset( bandNo ); + if ( myScale != 1.0 && myOffset != 0.0 ) + { + // if the band has scale and offset to apply change dataType + switch ( myDataType ) + { + case QGis::UnknownDataType: + case QGis::ARGB32: + case QGis::ARGB32_Premultiplied: + return myDataType; + break; + case QGis::Byte: + case QGis::UInt16: + case QGis::Int16: + case QGis::UInt32: + case QGis::Int32: + case QGis::Float32: + case QGis::CInt16: + myDataType = QGis::Float32; + break; + case QGis::Float64: + case QGis::CInt32: + case QGis::CFloat32: + myDataType = QGis::Float64; + break; + case QGis::CFloat64: + return myDataType; + break; + } + } + return myDataType; } QGis::DataType QgsGdalProvider::dataType( int bandNo ) const @@ -1125,6 +1117,38 @@ QGis::DataType QgsGdalProvider::dataType( int bandNo ) const return dataTypeFromGdal( mGdalDataType[bandNo-1] ); } +double QgsGdalProvider::bandScale( int bandNo ) const +{ +#if GDAL_VERSION_NUM >= 1800 + GDALRasterBandH myGdalBand = GDALGetRasterBand( mGdalDataset, bandNo ); + int bGotScale; + double myScale = GDALGetRasterScale( myGdalBand, &bGotScale ); + if ( bGotScale ) + return myScale; + else + return 1.0; +#else + Q_UNUSED( bandNo ); + return 1.0; +#endif +} + +double QgsGdalProvider::bandOffset( int bandNo ) const +{ +#if GDAL_VERSION_NUM >= 1800 + GDALRasterBandH myGdalBand = GDALGetRasterBand( mGdalDataset, bandNo ); + int bGotOffset; + double myOffset = GDALGetRasterOffset( myGdalBand, &bGotOffset ); + if ( bGotOffset ) + return myOffset; + else + return 0.0; +#else + Q_UNUSED( bandNo ); + return 0.0; +#endif +} + int QgsGdalProvider::bandCount() const { if ( mGdalDataset ) @@ -1355,10 +1379,19 @@ QgsRasterHistogram QgsGdalProvider::histogram( int theBandNo, // Min/max, if not specified, are set by histogramDefaults, it does not // set however min/max shifted to avoid rounding errors - + double myMinVal = myHistogram.minimum; double myMaxVal = myHistogram.maximum; + // unapply scale anf offset for min and max + double myScale = bandScale( theBandNo ); + double myOffset = bandOffset( theBandNo ); + if ( myScale != 1.0 || myOffset != 0. ) + { + myMinVal = (myHistogram.minimum - myOffset) / myScale; + myMaxVal = (myHistogram.maximum - myOffset) / myScale; + } + double dfHalfBucket = ( myMaxVal - myMinVal ) / ( 2 * myHistogram.binCount ); myMinVal -= dfHalfBucket; myMaxVal += dfHalfBucket; @@ -2352,6 +2385,35 @@ QgsRasterBandStats QgsGdalProvider::bandStatistics( int theBandNo, int theStats, myRasterBandStats.statsGathered = QgsRasterBandStats::Min | QgsRasterBandStats::Max | QgsRasterBandStats::Range | QgsRasterBandStats::Mean | QgsRasterBandStats::StdDev; + + // define if the band has scale and offset to apply + double myScale = bandScale( theBandNo ); + double myOffset = bandOffset( theBandNo ); + if ( myScale != 1.0 || myOffset != 0.0 ) + { + if ( myScale < 0.0 ) + { + // update Min and Max value + myRasterBandStats.minimumValue = pdfMax * myScale + myOffset; + myRasterBandStats.maximumValue = pdfMin * myScale + myOffset; + // update the range + myRasterBandStats.range = (pdfMin - pdfMax) * myScale; + // update standard deviation + myRasterBandStats.stdDev = -1.0 * pdfStdDev * myScale; + } + else + { + // update Min and Max value + myRasterBandStats.minimumValue = pdfMin * myScale + myOffset; + myRasterBandStats.maximumValue = pdfMax * myScale + myOffset; + // update the range + myRasterBandStats.range = (pdfMax - pdfMin) * myScale; + // update standard deviation + myRasterBandStats.stdDev = pdfStdDev * myScale; + } + // update the mean + myRasterBandStats.mean = pdfMean * myScale + myOffset; + } #ifdef QGISDEBUG QgsDebugMsg( "************ STATS **************" ); @@ -2562,6 +2624,36 @@ void QgsGdalProvider::initBaseDataset() #endif //mGdalDataType.append( myInternalGdalDataType ); + // define if the band has scale and offset to apply + double myScale = bandScale( i ); + double myOffset = bandOffset( i ); + if ( myScale != 1.0 && myOffset != 0.0 ) + { + // if the band has scale and offset to apply change dataType + switch ( myGdalDataType ) + { + case GDT_Unknown: + case GDT_TypeCount: + break; + case GDT_Byte: + case GDT_UInt16: + case GDT_Int16: + case GDT_UInt32: + case GDT_Int32: + case GDT_Float32: + case GDT_CInt16: + myGdalDataType = GDT_Float32; + break; + case GDT_Float64: + case GDT_CInt32: + case GDT_CFloat32: + myGdalDataType = GDT_Float64; + break; + case GDT_CFloat64: + break; + } + } + mGdalDataType.append( myGdalDataType ); //mInternalNoDataValue.append( myInternalNoDataValue ); //QgsDebugMsg( QString( "mInternalNoDataValue[%1] = %2" ).arg( i - 1 ).arg( mInternalNoDataValue[i-1] ) ); diff --git a/src/providers/gdal/qgsgdalprovider.h b/src/providers/gdal/qgsgdalprovider.h index ff9a22b42c22..30ebbdf4f304 100644 --- a/src/providers/gdal/qgsgdalprovider.h +++ b/src/providers/gdal/qgsgdalprovider.h @@ -157,8 +157,6 @@ class QgsGdalProvider : public QgsRasterDataProvider, QgsGdalProviderBase QGis::DataType dataType( int bandNo ) const; QGis::DataType srcDataType( int bandNo ) const; - QGis::DataType dataTypeFormGdal( int theGdalDataType ) const; - int bandCount() const; int colorInterpretation( int bandNo ) const; @@ -177,6 +175,13 @@ class QgsGdalProvider : public QgsRasterDataProvider, QgsGdalProviderBase void readBlock( int bandNo, int xBlock, int yBlock, void *data ); void readBlock( int bandNo, QgsRectangle const & viewExtent, int width, int height, void *data ); + /** Read band scale for raster value + * @@note added in 2.3 */ + double bandScale( int bandNo ) const; + /** Read band offset for raster value + * @@note added in 2.3 */ + double bandOffset( int bandNo ) const; + QList colorTable( int bandNo )const; /** diff --git a/tests/src/core/testqgsrasterlayer.cpp b/tests/src/core/testqgsrasterlayer.cpp index a03969831f33..429a2b43bf98 100644 --- a/tests/src/core/testqgsrasterlayer.cpp +++ b/tests/src/core/testqgsrasterlayer.cpp @@ -32,6 +32,7 @@ #include #include #include +#include #include #include #include @@ -66,6 +67,7 @@ class TestQgsRasterLayer: public QObject void landsatBasic875Qml(); void checkDimensions(); void checkStats(); + void checkScaleOffset(); void buildExternalOverviews(); void registry(); void transparency(); @@ -348,6 +350,87 @@ void TestQgsRasterLayer::checkStats() mReport += "

Passed

"; } +// test scale_factor and offset - uses netcdf file which may not be supported +// see http://hub.qgis.org/issues/8417 +void TestQgsRasterLayer::checkScaleOffset() +{ + mReport += "

Check Stats with scale/offset

\n"; + + QFileInfo myRasterFileInfo( mTestDataDir + "scaleoffset.tif" ); + QgsRasterLayer * myRasterLayer; + myRasterLayer = new QgsRasterLayer( myRasterFileInfo.filePath(), + myRasterFileInfo.completeBaseName() ); + QVERIFY ( myRasterLayer ); + if ( ! myRasterLayer->isValid() ) + { + qDebug() << QString( "raster layer %1 invalid" ).arg( myRasterFileInfo.filePath() ) ; + mReport += QString( "raster layer %1 invalid" ).arg( myRasterFileInfo.filePath() ) ; + delete myRasterLayer; + QVERIFY ( false ); + } + + QFile::remove( myRasterFileInfo.filePath() + ".aux.xml" ); // remove cached stats + QgsRasterBandStats myStatistics = myRasterLayer->dataProvider()->bandStatistics( 1, + QgsRasterBandStats::Min | QgsRasterBandStats::Max | + QgsRasterBandStats::Mean | QgsRasterBandStats::StdDev ); + mReport += QString( "raster min: %1 max: %2 mean: %3" ).arg( myStatistics.minimumValue ).arg( myStatistics.maximumValue ).arg( myStatistics.mean ); + QVERIFY( myRasterLayer->width() == 10 ); + QVERIFY( myRasterLayer->height() == 10 ); + //QVERIFY( myStatistics.elementCount == 100 ); + double minVal = 0.0; + mReport += QString( "min = %1 expected = %2 diff = %3
\n" ).arg( myStatistics.minimumValue ).arg( minVal ).arg( fabs( myStatistics.minimumValue - minVal ) ); + double maxVal = 9.0; + mReport += QString( "max = %1 expected = %2 diff = %3
\n" ).arg( myStatistics.maximumValue ).arg( maxVal ).arg( fabs( myStatistics.maximumValue - maxVal ) ); + double meanVal = 4.5; + mReport += QString( "min = %1 expected = %2 diff = %3
\n" ).arg( myStatistics.mean ).arg( meanVal ).arg( fabs( myStatistics.mean - meanVal ) ); + QVERIFY( fabs( myStatistics.minimumValue - minVal ) < 0.0000001 ); + QVERIFY( fabs( myStatistics.maximumValue - maxVal ) < 0.0000001 ); + QVERIFY( fabs( myStatistics.mean - meanVal ) < 0.0000001 ); + + double stdDev = 2.87228615; + // TODO: verify why GDAL stdDev is so different from generic (2.88675) + mReport += QString( "stdDev = %1 expected = %2 diff = %3
\n" ).arg( myStatistics.stdDev ).arg( stdDev ).arg( fabs( myStatistics.stdDev - stdDev ) ); + QVERIFY( fabs( myStatistics.stdDev - stdDev ) < 0.0000001 ); + + QgsRasterDataProvider* myProvider = myRasterLayer->dataProvider(); + QgsPoint myPoint(1535030, 5083350); + QgsRectangle myRect(1535030-5, 5083350-5, 1535030+5, 5083350+5); + QgsRasterIdentifyResult identifyResult = myProvider->identify( myPoint, QgsRaster::IdentifyFormatValue, myRect, 1, 1 ); + + if ( identifyResult.isValid() ) + { + QMap values = identifyResult.results(); + foreach ( int bandNo, values.keys() ) + { + QString valueString; + if ( values.value( bandNo ).isNull() ) + { + valueString = tr( "no data" ); + mReport += QString( " %1 = %2
\n" ).arg( myProvider->generateBandName( bandNo ) ).arg( valueString ); + delete myRasterLayer; + QVERIFY ( false ); + } + else + { + double expected = 0.99995432; + double value = values.value( bandNo ).toDouble(); + valueString = QgsRasterBlock::printValue( value ); + mReport += QString( " %1 = %2
\n" ).arg( myProvider->generateBandName( bandNo ) ).arg( valueString ); + mReport += QString( " value = %1 expected = %2 diff = %3
\n" ).arg( value ).arg( expected ).arg( fabs( value - expected ) ); + QVERIFY( fabs( value - expected ) < 0.0000001 ); + } + } + } + else + { + delete myRasterLayer; + QVERIFY ( false ); + } + + mReport += "

Passed

"; + delete myRasterLayer; +} + void TestQgsRasterLayer::buildExternalOverviews() { //before we begin delete any old ovr file (if it exists) diff --git a/tests/testdata/scaleoffset.tif b/tests/testdata/scaleoffset.tif new file mode 100644 index 000000000000..713a52f9ccd6 Binary files /dev/null and b/tests/testdata/scaleoffset.tif differ