From 07c57585c1fe10143a3529dfdeddcfa358c8f692 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=27Hont=20Ren=C3=A9-Luc?= Date: Fri, 10 Jan 2014 10:36:12 +0100 Subject: [PATCH] [RASTER][Feature] Applying scale and offset to raster data - funded Ifremer An issue has been opened 5 mounth ago [BUG] #8417 incorrect value loaded from netcdf file The data has not be loaded incorrectly, but QGIS doesn't apply scale and offset defined for each band. This commit will apply scale and offset to GDAL Provider BandStatistics and to RASTER block of data. It also adds bandScale and bandOffset method to QgsRasterDataProvider Python API. --- python/core/raster/qgsrasterblock.sip | 4 + python/core/raster/qgsrasterdataprovider.sip | 8 + src/core/raster/qgsrasterblock.cpp | 14 ++ src/core/raster/qgsrasterblock.h | 4 + src/core/raster/qgsrasterdataprovider.cpp | 2 + src/core/raster/qgsrasterdataprovider.h | 7 + src/providers/gdal/qgsgdalprovider.cpp | 184 ++++++++++++++----- src/providers/gdal/qgsgdalprovider.h | 9 +- tests/src/core/testqgsrasterlayer.cpp | 83 +++++++++ tests/testdata/scaleoffset.tif | Bin 0 -> 2042 bytes 10 files changed, 267 insertions(+), 48 deletions(-) create mode 100644 tests/testdata/scaleoffset.tif 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 0000000000000000000000000000000000000000..713a52f9ccd669779f1884cc15437357c96824c5 GIT binary patch literal 2042 zcmdT_%Wl&^6rH4~psFgc02Z(qTqN46ar{csRH=pQB$ZJfs5TKhYlx?9gfq762~t;- z$Cf`}$%-#viP*po@CUGA#}6RJt|3kvC*5KxnKRe-+RPLfDRK}+>Jo%CgW@41Xoz>XT|^vl z)ou=4w};#j=F~!!J+)%wGu;_Ej_n;BgqU)w!k1>bx*j>SqePqIF=azzc#P$`x6zC|lEcE~2ZcbZk1HBL`NY87Kc?rBvu=}OZSqjN4<;bPfLt6p#PCh|!_RaAJP zB9%&8RVm|!>^7Tb!%FjD;6)`c_ii9t~10OJn= z0P&T7`mq4_9-M-|V`9!2?z