[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.
qgis · May 15, 2014 · 07c5758 · 07c5758
1 parent 60f7824
commit 07c5758
Show file tree

Hide file tree

Showing 10 changed files with 267 additions and 48 deletions.
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;
 
     /**