sharpetronics/QGIS
1
0
Fork 0
mirror of https://github.com/qgis/QGIS.git synced 2025-12-15 00:07:25 -05:00
Code Issues Packages Projects Releases Wiki Activity

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

Browse Source 
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.
This commit is contained in:
D'Hont René-Luc 2014-01-10 10:36:12 +01:00
parent 60f782441a
commit 07c57585c1
10 changed files with 267 additions and 48 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
python/core/raster/qgsrasterblock.sip
View File

@ -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;

8
python/core/raster/qgsrasterdataprovider.sip
View File

@ -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;

14
src/core/raster/qgsrasterblock.cpp
View File

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

4
src/core/raster/qgsrasterblock.h
View File

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

2
src/core/raster/qgsrasterdataprovider.cpp
View File

@ -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;

7
src/core/raster/qgsrasterdataprovider.h
View File

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

184
src/providers/gdal/qgsgdalprovider.cpp
View File

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

9
src/providers/gdal/qgsgdalprovider.h
View File

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

83
tests/src/core/testqgsrasterlayer.cpp
View File

@ -32,6 +32,7 @@
#include <qgsrasterpyramid.h>
#include <qgsrasterbandstats.h>
#include <qgsrasterpyramid.h>
#include <qgsrasteridentifyresult.h>
#include <qgsmaplayerregistry.h>
#include <qgsapplication.h>
#include <qgsmaprenderer.h>
@ -66,6 +67,7 @@ class TestQgsRasterLayer: public QObject
void landsatBasic875Qml();
void checkDimensions();
void checkStats();
void checkScaleOffset();
void buildExternalOverviews();
void registry();
void transparency();
@ -348,6 +350,87 @@ 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::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)

BIN
tests/testdata/scaleoffset.tif vendored Normal file
View File

Binary file not shown.