sharpetronics/QGIS
1
0
Fork 0
mirror of https://github.com/qgis/QGIS.git synced 2025-10-16 00:05:45 -04:00
Code Issues Packages Projects Releases Wiki Activity

[rastercalc] More robust handling of nodata in calculations

Browse Source 
Also allow creation of QgsRasterCalcNodes which directly reference
a QgsRasterMatrix for testing.
This commit is contained in:
Nyall Dawson 2015-06-04 20:29:02 +10:00
parent 1219a0fbd2
commit e1f7d330c2
4 changed files with 251 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
python/analysis/raster/qgsrastercalcnode.sip
View File

@ -9,7 +9,8 @@ class QgsRasterCalcNode
{
tOperator,
tNumber,
tRasterRef
tRasterRef,
tMatrix
};
//! possible operators
@ -41,6 +42,7 @@ class QgsRasterCalcNode
QgsRasterCalcNode();
QgsRasterCalcNode( double number );
QgsRasterCalcNode( QgsRasterMatrix* matrix );
QgsRasterCalcNode( Operator op, QgsRasterCalcNode* left, QgsRasterCalcNode* right );
QgsRasterCalcNode( const QString& rasterName );
~QgsRasterCalcNode();

34
src/analysis/raster/qgsrastercalcnode.cpp
View File

@ -20,6 +20,7 @@ QgsRasterCalcNode::QgsRasterCalcNode()
, mLeft( 0 )
, mRight( 0 )
, mNumber( 0 )
, mMatrix( 0 )
, mOperator( opNONE )
{
}
@ -29,15 +30,28 @@ QgsRasterCalcNode::QgsRasterCalcNode( double number )
, mLeft( 0 )
, mRight( 0 )
, mNumber( number )
, mMatrix( 0 )
, mOperator( opNONE )
{
}
QgsRasterCalcNode::QgsRasterCalcNode( QgsRasterMatrix* matrix )
: mType( tMatrix )
, mLeft( 0 )
, mRight( 0 )
, mNumber( 0 )
, mMatrix( matrix )
, mOperator( opNONE )
{
}
QgsRasterCalcNode::QgsRasterCalcNode( Operator op, QgsRasterCalcNode* left, QgsRasterCalcNode* right )
: mType( tOperator )
, mLeft( left )
, mRight( right )
, mNumber( 0 )
, mMatrix( 0 )
, mOperator( op )
{
}
@ -48,6 +62,7 @@ QgsRasterCalcNode::QgsRasterCalcNode( const QString& rasterName )
, mRight( 0 )
, mNumber( 0 )
, mRasterName( rasterName )
, mMatrix( 0 )
, mOperator( opNONE )
{
if ( mRasterName.startsWith( '"' ) && mRasterName.endsWith( '"' ) )
@ -81,8 +96,12 @@ bool QgsRasterCalcNode::calculate( QMap<QString, QgsRasterMatrix*>& rasterData,
int nEntries = ( *it )->nColumns() * ( *it )->nRows();
double* data = new double[nEntries];
memcpy( data, ( *it )->data(), nEntries * sizeof( double ) );
result.setData(( *it )->nColumns(), ( *it )->nRows(), data, ( *it )->nodataValue() );
for ( int i = 0; i < nEntries; ++i )
{
data[i] = ( *it )->data()[i] == ( *it )->nodataValue() ? result.nodataValue() : ( *it )->data()[i];
}
result.setData(( *it )->nColumns(), ( *it )->nRows(), data, result.nodataValue() );
return true;
}
else if ( mType == tOperator )
@ -186,6 +205,17 @@ bool QgsRasterCalcNode::calculate( QMap<QString, QgsRasterMatrix*>& rasterData,
result.setData( 1, 1, data, result.nodataValue() );
return true;
}
else if ( mType == tMatrix )
{
int nEntries = mMatrix->nColumns() * mMatrix->nRows();
double* data = new double[nEntries];
for ( int i = 0; i < nEntries; ++i )
{
data[i] = mMatrix->data()[i] == mMatrix->nodataValue() ? result.nodataValue() : mMatrix->data()[i];
}
result.setData( mMatrix->nColumns(), mMatrix->nRows(), data, result.nodataValue() );
return true;
}
return false;
}

5
src/analysis/raster/qgsrastercalcnode.h
View File

@ -31,7 +31,8 @@ class ANALYSIS_EXPORT QgsRasterCalcNode
{
tOperator = 1,
tNumber,
tRasterRef
tRasterRef,
tMatrix
};
//! possible operators
@ -65,6 +66,7 @@ class ANALYSIS_EXPORT QgsRasterCalcNode
QgsRasterCalcNode();
QgsRasterCalcNode( double number );
QgsRasterCalcNode( QgsRasterMatrix* matrix );
QgsRasterCalcNode( Operator op, QgsRasterCalcNode* left, QgsRasterCalcNode* right );
QgsRasterCalcNode( const QString& rasterName );
~QgsRasterCalcNode();
@ -86,6 +88,7 @@ class ANALYSIS_EXPORT QgsRasterCalcNode
QgsRasterCalcNode* mRight;
double mNumber;
QString mRasterName;
QgsRasterMatrix* mMatrix;
Operator mOperator;
};

212
tests/src/analysis/testqgsrastercalculator.cpp
View File

@ -40,6 +40,14 @@ class TestQgsRasterCalculator : public QObject
void singleOp_data();
void singleOp(); //test operators which operate on a single value
void singleOpMatrices(); // test single op using matrix
void dualOpNumberMatrix(); // test dual op run on number and matrix
void dualOpMatrixNumber(); // test dual op run on matrix and number
void dualOpMatrixMatrix(); // test dual op run on matrix and matrix
void rasterRefOp();
void dualOpRasterRaster(); //test dual op on raster ref and raster ref
};
void TestQgsRasterCalculator::initTestCase()
@ -178,5 +186,209 @@ void TestQgsRasterCalculator::singleOp()
}
void TestQgsRasterCalculator::singleOpMatrices()
{
// test single op run on matrix
double* d = new double[6];
d[0] = 1.0;
d[1] = 2.0;
d[2] = 3.0;
d[3] = 4.0;
d[4] = 5.0;
d[5] = -1.0;
QgsRasterMatrix m( 2, 3, d, -1.0 );
QgsRasterCalcNode node( QgsRasterCalcNode::opSIGN, new QgsRasterCalcNode( &m ) , 0 );
QgsRasterMatrix result;
result.setNodataValue( -9999 );
QMap<QString, QgsRasterMatrix*> rasterData;
QVERIFY( node.calculate( rasterData, result ) );
QCOMPARE( result.data()[0], -d[0] );
QCOMPARE( result.data()[1], -d[1] );
QCOMPARE( result.data()[2], -d[2] );
QCOMPARE( result.data()[3], -d[3] );
QCOMPARE( result.data()[4], -d[4] );
QCOMPARE( result.data()[5], -9999.0 );
}
void TestQgsRasterCalculator::dualOpNumberMatrix()
{
// test dual op run on number and matrix
double* d = new double[6];
d[0] = 1.0;
d[1] = 2.0;
d[2] = 3.0;
d[3] = 4.0;
d[4] = 5.0;
d[5] = -1.0;
QgsRasterMatrix m( 2, 3, d, -1.0 );
QgsRasterCalcNode node( QgsRasterCalcNode::opPLUS, new QgsRasterCalcNode( 5.0 ), new QgsRasterCalcNode( &m ) );
QgsRasterMatrix result;
result.setNodataValue( -9999 );
QMap<QString, QgsRasterMatrix*> rasterData;
QVERIFY( node.calculate( rasterData, result ) );
QCOMPARE( result.data()[0], 6.0 );
QCOMPARE( result.data()[1], 7.0 );
QCOMPARE( result.data()[2], 8.0 );
QCOMPARE( result.data()[3], 9.0 );
QCOMPARE( result.data()[4], 10.0 );
QCOMPARE( result.data()[5], -9999.0 );
//also check adding no data number
QgsRasterCalcNode nodeNoData( QgsRasterCalcNode::opPLUS, new QgsRasterCalcNode( -9999 ), new QgsRasterCalcNode( &m ) );
QVERIFY( nodeNoData.calculate( rasterData, result ) );
QCOMPARE( result.data()[0], -9999.0 );
}
void TestQgsRasterCalculator::dualOpMatrixNumber()
{
// test dual op run on matrix and number
double* d = new double[6];
d[0] = 1.0;
d[1] = 2.0;
d[2] = 3.0;
d[3] = 4.0;
d[4] = 5.0;
d[5] = -1.0;
QgsRasterMatrix m( 2, 3, d, -1.0 );
QgsRasterCalcNode node( QgsRasterCalcNode::opPLUS, new QgsRasterCalcNode( &m ), new QgsRasterCalcNode( 5.0 ) );
QgsRasterMatrix result;
result.setNodataValue( -9999 );
QMap<QString, QgsRasterMatrix*> rasterData;
QVERIFY( node.calculate( rasterData, result ) );
QCOMPARE( result.data()[0], 6.0 );
QCOMPARE( result.data()[1], 7.0 );
QCOMPARE( result.data()[2], 8.0 );
QCOMPARE( result.data()[3], 9.0 );
QCOMPARE( result.data()[4], 10.0 );
QCOMPARE( result.data()[5], -9999.0 );
//also check adding no data number
QgsRasterCalcNode nodeNoData( QgsRasterCalcNode::opPLUS, new QgsRasterCalcNode( &m ), new QgsRasterCalcNode( -9999 ) );
QVERIFY( nodeNoData.calculate( rasterData, result ) );
QCOMPARE( result.data()[0], -9999.0 );
}
void TestQgsRasterCalculator::dualOpMatrixMatrix()
{
// test dual op run on matrix and matrix
double* d = new double[6];
d[0] = 1.0;
d[1] = 2.0;
d[2] = -2.0;
d[3] = -1.0; //nodata
d[4] = 5.0;
d[5] = -1.0; //nodata
QgsRasterMatrix m1( 2, 3, d, -1.0 );
double* d2 = new double[6];
d2[0] = -1.0;
d2[1] = -2.0; //nodata
d2[2] = 13.0;
d2[3] = -2.0; //nodata
d2[4] = 15.0;
d2[5] = -1.0;
QgsRasterMatrix m2( 2, 3, d2, -2.0 ); //different no data value
QgsRasterCalcNode node( QgsRasterCalcNode::opPLUS, new QgsRasterCalcNode( &m1 ), new QgsRasterCalcNode( &m2 ) );
QgsRasterMatrix result;
result.setNodataValue( -9999 );
QMap<QString, QgsRasterMatrix*> rasterData;
QVERIFY( node.calculate( rasterData, result ) );
QCOMPARE( result.data()[0], 0.0 );
QCOMPARE( result.data()[1], -9999.0 );
QCOMPARE( result.data()[2], 11.0 );
QCOMPARE( result.data()[3], -9999.0 );
QCOMPARE( result.data()[4], 20.0 );
QCOMPARE( result.data()[5], -9999.0 );
}
void TestQgsRasterCalculator::rasterRefOp()
{
// test single op run on raster ref
QgsRasterCalcNode node( QgsRasterCalcNode::opSIGN, new QgsRasterCalcNode( "raster" ), 0 );
QgsRasterMatrix result;
result.setNodataValue( -9999 );
QMap<QString, QgsRasterMatrix*> rasterData;
//first test invalid raster ref
QVERIFY( !node.calculate( rasterData, result ) );
//now create raster ref
double* d = new double[6];
d[0] = 1.0;
d[1] = 2.0;
d[2] = 3.0;
d[3] = 4.0;
d[4] = 5.0;
d[5] = -1.0;
QgsRasterMatrix m( 2, 3, d, -1.0 );
rasterData.insert( "raster", &m );
QVERIFY( node.calculate( rasterData, result ) );
QCOMPARE( result.data()[0], -d[0] );
QCOMPARE( result.data()[1], -d[1] );
QCOMPARE( result.data()[2], -d[2] );
QCOMPARE( result.data()[3], -d[3] );
QCOMPARE( result.data()[4], -d[4] );
QCOMPARE( result.data()[5], -9999.0 );
}
void TestQgsRasterCalculator::dualOpRasterRaster()
{
// test dual op run on matrix and matrix
double* d = new double[6];
d[0] = 1.0;
d[1] = 2.0;
d[2] = -2.0;
d[3] = -1.0; //nodata
d[4] = 5.0;
d[5] = -1.0; //nodata
QgsRasterMatrix m1( 2, 3, d, -1.0 );
QMap<QString, QgsRasterMatrix*> rasterData;
rasterData.insert( "raster1", &m1 );
double* d2 = new double[6];
d2[0] = -1.0;
d2[1] = -2.0; //nodata
d2[2] = 13.0;
d2[3] = -2.0; //nodata
d2[4] = 15.0;
d2[5] = -1.0;
QgsRasterMatrix m2( 2, 3, d2, -2.0 ); //different no data value
rasterData.insert( "raster2", &m2 );
QgsRasterCalcNode node( QgsRasterCalcNode::opPLUS, new QgsRasterCalcNode( "raster1" ), new QgsRasterCalcNode( "raster2" ) );
QgsRasterMatrix result;
result.setNodataValue( -9999 );
QVERIFY( node.calculate( rasterData, result ) );
QCOMPARE( result.data()[0], 0.0 );
QCOMPARE( result.data()[1], -9999.0 );
QCOMPARE( result.data()[2], 11.0 );
QCOMPARE( result.data()[3], -9999.0 );
QCOMPARE( result.data()[4], 20.0 );
QCOMPARE( result.data()[5], -9999.0 );
}
QTEST_MAIN( TestQgsRasterCalculator )
#include "testqgsrastercalculator.moc"