Skip to content
Permalink
Browse files
[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.
  • Loading branch information
rldhont committed May 15, 2014
1 parent 60f7824 commit 07c57585c1fe10143a3529dfdeddcfa358c8f692
@@ -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;

@@ -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;
@@ -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() )
@@ -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; }

@@ -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;
@@ -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. */
@@ -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] ) );
@@ -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;

/**

0 comments on commit 07c5758

Please sign in to comment.