Skip to content

Commit

Permalink
Merge branch 'raster_scale_offset'
Browse files Browse the repository at this point in the history
  • Loading branch information
@rldhont
rldhont committed Feb 6, 2014
2 parents d26831b + 77d0a7c commit bf3fcec
Show file tree
Hide file tree
Showing 8 changed files with 197 additions and 48 deletions.
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.1 */
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 @@ -367,6 +367,13 @@ class CORE_EXPORT QgsRasterDataProvider : public QgsDataProvider, public QgsRast
virtual void readBlock( int bandNo, QgsRectangle const & viewExtent, int width, int height, void *data )
{ Q_UNUSED( bandNo ); Q_UNUSED( viewExtent ); Q_UNUSED( width ); Q_UNUSED( height ); Q_UNUSED( data ); }

/** Read band scale for raster value
* @@note added in 2.1 */
virtual double bandScale( int bandNo ) { Q_UNUSED( bandNo ); return 1.0; }
/** Read band offset for raster value
* @@note added in 2.1 */
virtual double bandOffset( int bandNo ) { Q_UNUSED( bandNo ); return 0.0; }

/** Returns true if user no data contains value */
bool userNoDataValuesContains( int bandNo, double value ) const;

Expand Down
168 changes: 122 additions & 46 deletions src/providers/gdal/qgsgdalprovider.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -386,6 +386,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 @@ -994,55 +996,43 @@ 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:
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 @@ -1052,6 +1042,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 @@ -1282,10 +1304,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 @@ -2252,6 +2283,22 @@ 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 )
{
// update Min and MAx value
myRasterBandStats.minimumValue = pdfMin * myScale + myOffset;
myRasterBandStats.maximumValue = pdfMax * myScale + myOffset;
// update the range
myRasterBandStats.range = (pdfMax - pdfMin) * myScale;
// update the mean
myRasterBandStats.mean = pdfMean * myScale + myOffset;
// update standard deviation
myRasterBandStats.stdDev = pdfStdDev * myScale;
}

#ifdef QGISDEBUG
QgsDebugMsg( "************ STATS **************" );
Expand Down Expand Up @@ -2463,6 +2510,35 @@ 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:
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 @@ -175,6 +173,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.1 */
double bandScale( int bandNo ) const;
/** Read band offset for raster value
* @@note added in 2.1 */
double bandOffset( int bandNo ) const;

QList<QgsColorRampShader::ColorRampItem> colorTable( int bandNo )const;

/**
Expand Down
41 changes: 41 additions & 0 deletions tests/src/core/testqgsrasterlayer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,7 @@ class TestQgsRasterLayer: public QObject
void landsatBasic875Qml();
void checkDimensions();
void checkStats();
void checkStatsScaleOffset();
void buildExternalOverviews();
void registry();
void transparency();
Expand Down Expand Up @@ -355,6 +356,46 @@ void TestQgsRasterLayer::checkStats()
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::checkStatsScaleOffset()
{
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 mpRasterLayer;
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 );
qDebug() << 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 );
QVERIFY( myStatistics.minimumValue == 0 );
QVERIFY( myStatistics.maximumValue == 9 );
QVERIFY( myStatistics.mean == 4.5 );
double stdDev = 2.87228132326901431;
// TODO: verify why GDAL stdDev is so different from generic (2.88675)
mReport += QString( "stdDev = %1 expected = %2<br>\n" ).arg( myStatistics.stdDev ).arg( stdDev );
QVERIFY( fabs( myStatistics.stdDev - stdDev )
< 0.0000000000000001 );
mReport += "<p>Passed</p>";

delete myRasterLayer;
}

void TestQgsRasterLayer::buildExternalOverviews()
{
//before we begin delete any old ovr file (if it exists)
Expand Down
Binary file added tests/testdata/scaleoffset.tif
View file
Binary file not shown.

0 comments on commit bf3fcec

Please sign in to comment.