Skip to content

Commit

Permalink
[rastercalc] More robust handling of nodata in calculations
Browse files Browse the repository at this point in the history 
Also allow creation of QgsRasterCalcNodes which directly reference
a QgsRasterMatrix for testing.
  • Loading branch information
@nyalldawson
nyalldawson committed Jun 10, 2015
1 parent 1219a0f commit e1f7d33
Show file tree
Hide file tree
Showing 4 changed files with 251 additions and 4 deletions.
4 changes: 3 additions & 1 deletion python/analysis/raster/qgsrastercalcnode.sip
View file
Expand Up @@ -9,7 +9,8 @@ class QgsRasterCalcNode
{
tOperator,
tNumber,
tRasterRef
tRasterRef,
tMatrix
};

//! possible operators
Expand Down Expand Up @@ -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();
Expand Down
34 changes: 32 additions & 2 deletions src/analysis/raster/qgsrastercalcnode.cpp
View file
Expand Up @@ -20,6 +20,7 @@ QgsRasterCalcNode::QgsRasterCalcNode()
, mLeft( 0 )
, mRight( 0 )
, mNumber( 0 )
, mMatrix( 0 )
, mOperator( opNONE )
{
}
Expand All @@ -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 )
{
}
Expand All @@ -48,6 +62,7 @@ QgsRasterCalcNode::QgsRasterCalcNode( const QString& rasterName )
, mRight( 0 )
, mNumber( 0 )
, mRasterName( rasterName )
, mMatrix( 0 )
, mOperator( opNONE )
{
if ( mRasterName.startsWith( '"' ) && mRasterName.endsWith( '"' ) )
Expand Down Expand Up @@ -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 )
Expand Down Expand Up @@ -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;
}

Expand Down
5 changes: 4 additions & 1 deletion src/analysis/raster/qgsrastercalcnode.h
View file
Expand Up @@ -31,7 +31,8 @@ class ANALYSIS_EXPORT QgsRasterCalcNode
{
tOperator = 1,
tNumber,
tRasterRef
tRasterRef,
tMatrix
};

//! possible operators
Expand Down Expand Up @@ -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();
Expand All @@ -86,6 +88,7 @@ class ANALYSIS_EXPORT QgsRasterCalcNode
QgsRasterCalcNode* mRight;
double mNumber;
QString mRasterName;
QgsRasterMatrix* mMatrix;
Operator mOperator;
};

Expand Down
212 changes: 212 additions & 0 deletions tests/src/analysis/testqgsrastercalculator.cpp
View file
Expand Up @@ -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()
Expand Down Expand Up @@ -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"

0 comments on commit e1f7d33

Please sign in to comment.