Skip to content

Commit

Permalink
[RASTER][Feature] Applying scale and offset to raster data - funded
Browse files Browse the repository at this point in the history 
Ifremer

An issue has been opened 5 mounth ago
[BUG] qgis#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
rldhont committed May 15, 2014
1 parent 60f7824 commit 07c5758
Show file tree
Hide file tree
Showing 10 changed files with 267 additions and 48 deletions.
4 changes: 4 additions & 0 deletions python/core/raster/qgsrasterblock.sip
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down
8 changes: 8 additions & 0 deletions python/core/raster/qgsrasterdataprovider.sip
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
14 changes: 14 additions & 0 deletions src/core/raster/qgsrasterblock.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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() )
Expand Down
4 changes: 4 additions & 0 deletions src/core/raster/qgsrasterblock.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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; }

Expand Down
2 changes: 2 additions & 0 deletions src/core/raster/qgsrasterdataprovider.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
7 changes: 7 additions & 0 deletions src/core/raster/qgsrasterdataprovider.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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. */
Expand Down
184 changes: 138 additions & 46 deletions src/providers/gdal/qgsgdalprovider.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
Expand Down Expand Up @@ -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
Expand All @@ -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 )
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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 **************" );
Expand Down Expand Up @@ -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] ) );
Expand Down
9 changes: 7 additions & 2 deletions src/providers/gdal/qgsgdalprovider.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;

/**
Expand Down
Loading

0 comments on commit 07c5758

Please sign in to comment.