[RASTER][Feature] Applying scale and offset to raster data - funded

rldhont · rldhont · commit 07c57585c1fe · 2014-05-15T16:55:30.000+02:00
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.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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<QgsColorRampShader::ColorRampItem> colorTable( int bandNo )const;
 
     /**
diff --git a/tests/src/core/testqgsrasterlayer.cpp b/tests/src/core/testqgsrasterlayer.cpp
diff --git a/tests/testdata/scaleoffset.tif b/tests/testdata/scaleoffset.tif

Original file line number	Diff line number	Diff line change
`@@ -208,6 +208,8 @@ QgsRasterBlock * QgsRasterDataProvider::block( int theBandNo, QgsRectangle cons`
`208`	`208`	`readBlock( theBandNo, theExtent, theWidth, theHeight, block->bits() );`
`209`	`209`	`}`
`210`	`210`
	`211`	`+ // apply scale and offset`
	`212`	`+ block->applyScaleOffset( bandScale( theBandNo ), bandOffset( theBandNo ) );`
`211`	`213`	`// apply user no data values`
`212`	`214`	`block->applyNoDataValues( userNoDataValues( theBandNo ) );`
`213`	`215`	`return block;`
-Original file line number
+Diff line change
 #include <qgsrasterpyramid.h>
 #include <qgsrasterbandstats.h>
 #include <qgsrasterpyramid.h>
 +#include <qgsrasteridentifyresult.h>
 #include <qgsmaplayerregistry.h>
 #include <qgsapplication.h>
 #include <qgsmaprenderer.h>
     void landsatBasic875Qml();
     void checkDimensions();
     void checkStats();
 +    void checkScaleOffset();
     void buildExternalOverviews();
     void registry();
     void transparency();
   mReport += "<p>Passed</p>";
+}
 +// test scale_factor and offset - uses netcdf file which may not be supported
 +// see http://hub.qgis.org/issues/8417
 +void TestQgsRasterLayer::checkScaleOffset()
 +{
 +  mReport += "<h2>Check Stats with scale/offset</h2>\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<br>\n" ).arg( myStatistics.minimumValue ).arg( minVal ).arg( fabs( myStatistics.minimumValue - minVal ) );
 +  double maxVal = 9.0;
 +  mReport += QString( "max = %1 expected = %2 diff = %3<br>\n" ).arg( myStatistics.maximumValue ).arg( maxVal ).arg( fabs( myStatistics.maximumValue - maxVal ) );
 +  double meanVal = 4.5;
 +  mReport += QString( "min = %1 expected = %2 diff = %3<br>\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<br>\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<int, QVariant> values = identifyResult.results();
 +    foreach ( int bandNo, values.keys() )
 +    {
 +      QString valueString;
 +      if ( values.value( bandNo ).isNull() )
 +      {
 +        valueString = tr( "no data" );
 +        mReport += QString( " %1 = %2 <br>\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 <br>\n" ).arg( myProvider->generateBandName( bandNo ) ).arg( valueString );
 +        mReport += QString( " value = %1 expected = %2 diff = %3 <br>\n" ).arg( value ).arg( expected ).arg( fabs( value - expected ) );
 +        QVERIFY( fabs( value - expected ) < 0.0000001 );
 +      }
 +    }
 +  }
 +  else
 +  {
 +    delete myRasterLayer;
 +    QVERIFY ( false );
 +  }
++
 +  mReport += "<p>Passed</p>";
 +  delete myRasterLayer;
 +}
++
 void TestQgsRasterLayer::buildExternalOverviews()
+{
   //before we begin delete any old ovr file (if it exists)