Skip to content

Commit e1f7d33

Browse files
nyalldawsonnyalldawson
committedJun 10, 2015
[rastercalc] More robust handling of nodata in calculations
Also allow creation of QgsRasterCalcNodes which directly reference a QgsRasterMatrix for testing.
·
ltr-3_40archive/master_2
1 parent 1219a0f commit e1f7d33

File tree

4 files changed

+251
-4
lines changed

4 files changed

+251
-4
lines changed
 

‎python/analysis/raster/qgsrastercalcnode.sip

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,8 @@ class QgsRasterCalcNode
99
{
1010
tOperator,
1111
tNumber,
12-
tRasterRef
12+
tRasterRef,
13+
tMatrix
1314
};
1415

1516
//! possible operators
@@ -41,6 +42,7 @@ class QgsRasterCalcNode
4142

4243
QgsRasterCalcNode();
4344
QgsRasterCalcNode( double number );
45+
QgsRasterCalcNode( QgsRasterMatrix* matrix );
4446
QgsRasterCalcNode( Operator op, QgsRasterCalcNode* left, QgsRasterCalcNode* right );
4547
QgsRasterCalcNode( const QString& rasterName );
4648
~QgsRasterCalcNode();

‎src/analysis/raster/qgsrastercalcnode.cpp

Lines changed: 32 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@ QgsRasterCalcNode::QgsRasterCalcNode()
2020
, mLeft( 0 )
2121
, mRight( 0 )
2222
, mNumber( 0 )
23+
, mMatrix( 0 )
2324
, mOperator( opNONE )
2425
{
2526
}
@@ -29,15 +30,28 @@ QgsRasterCalcNode::QgsRasterCalcNode( double number )
2930
, mLeft( 0 )
3031
, mRight( 0 )
3132
, mNumber( number )
33+
, mMatrix( 0 )
3234
, mOperator( opNONE )
3335
{
3436
}
3537

38+
QgsRasterCalcNode::QgsRasterCalcNode( QgsRasterMatrix* matrix )
39+
: mType( tMatrix )
40+
, mLeft( 0 )
41+
, mRight( 0 )
42+
, mNumber( 0 )
43+
, mMatrix( matrix )
44+
, mOperator( opNONE )
45+
{
46+
47+
}
48+
3649
QgsRasterCalcNode::QgsRasterCalcNode( Operator op, QgsRasterCalcNode* left, QgsRasterCalcNode* right )
3750
: mType( tOperator )
3851
, mLeft( left )
3952
, mRight( right )
4053
, mNumber( 0 )
54+
, mMatrix( 0 )
4155
, mOperator( op )
4256
{
4357
}
@@ -48,6 +62,7 @@ QgsRasterCalcNode::QgsRasterCalcNode( const QString& rasterName )
4862
, mRight( 0 )
4963
, mNumber( 0 )
5064
, mRasterName( rasterName )
65+
, mMatrix( 0 )
5166
, mOperator( opNONE )
5267
{
5368
if ( mRasterName.startsWith( '"' ) && mRasterName.endsWith( '"' ) )
@@ -81,8 +96,12 @@ bool QgsRasterCalcNode::calculate( QMap<QString, QgsRasterMatrix*>& rasterData,
8196

8297
int nEntries = ( *it )->nColumns() * ( *it )->nRows();
8398
double* data = new double[nEntries];
84-
memcpy( data, ( *it )->data(), nEntries * sizeof( double ) );
85-
result.setData(( *it )->nColumns(), ( *it )->nRows(), data, ( *it )->nodataValue() );
99+
100+
for ( int i = 0; i < nEntries; ++i )
101+
{
102+
data[i] = ( *it )->data()[i] == ( *it )->nodataValue() ? result.nodataValue() : ( *it )->data()[i];
103+
}
104+
result.setData(( *it )->nColumns(), ( *it )->nRows(), data, result.nodataValue() );
86105
return true;
87106
}
88107
else if ( mType == tOperator )
@@ -186,6 +205,17 @@ bool QgsRasterCalcNode::calculate( QMap<QString, QgsRasterMatrix*>& rasterData,
186205
result.setData( 1, 1, data, result.nodataValue() );
187206
return true;
188207
}
208+
else if ( mType == tMatrix )
209+
{
210+
int nEntries = mMatrix->nColumns() * mMatrix->nRows();
211+
double* data = new double[nEntries];
212+
for ( int i = 0; i < nEntries; ++i )
213+
{
214+
data[i] = mMatrix->data()[i] == mMatrix->nodataValue() ? result.nodataValue() : mMatrix->data()[i];
215+
}
216+
result.setData( mMatrix->nColumns(), mMatrix->nRows(), data, result.nodataValue() );
217+
return true;
218+
}
189219
return false;
190220
}
191221

‎src/analysis/raster/qgsrastercalcnode.h

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -31,7 +31,8 @@ class ANALYSIS_EXPORT QgsRasterCalcNode
3131
{
3232
tOperator = 1,
3333
tNumber,
34-
tRasterRef
34+
tRasterRef,
35+
tMatrix
3536
};
3637

3738
//! possible operators
@@ -65,6 +66,7 @@ class ANALYSIS_EXPORT QgsRasterCalcNode
6566

6667
QgsRasterCalcNode();
6768
QgsRasterCalcNode( double number );
69+
QgsRasterCalcNode( QgsRasterMatrix* matrix );
6870
QgsRasterCalcNode( Operator op, QgsRasterCalcNode* left, QgsRasterCalcNode* right );
6971
QgsRasterCalcNode( const QString& rasterName );
7072
~QgsRasterCalcNode();
@@ -86,6 +88,7 @@ class ANALYSIS_EXPORT QgsRasterCalcNode
8688
QgsRasterCalcNode* mRight;
8789
double mNumber;
8890
QString mRasterName;
91+
QgsRasterMatrix* mMatrix;
8992
Operator mOperator;
9093
};
9194

‎tests/src/analysis/testqgsrastercalculator.cpp

Lines changed: 212 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,14 @@ class TestQgsRasterCalculator : public QObject
4040
void singleOp_data();
4141
void singleOp(); //test operators which operate on a single value
4242

43+
void singleOpMatrices(); // test single op using matrix
44+
void dualOpNumberMatrix(); // test dual op run on number and matrix
45+
void dualOpMatrixNumber(); // test dual op run on matrix and number
46+
void dualOpMatrixMatrix(); // test dual op run on matrix and matrix
47+
48+
void rasterRefOp();
49+
void dualOpRasterRaster(); //test dual op on raster ref and raster ref
50+
4351
};
4452

4553
void TestQgsRasterCalculator::initTestCase()
@@ -178,5 +186,209 @@ void TestQgsRasterCalculator::singleOp()
178186

179187
}
180188

189+
void TestQgsRasterCalculator::singleOpMatrices()
190+
{
191+
// test single op run on matrix
192+
double* d = new double[6];
193+
d[0] = 1.0;
194+
d[1] = 2.0;
195+
d[2] = 3.0;
196+
d[3] = 4.0;
197+
d[4] = 5.0;
198+
d[5] = -1.0;
199+
200+
QgsRasterMatrix m( 2, 3, d, -1.0 );
201+
202+
QgsRasterCalcNode node( QgsRasterCalcNode::opSIGN, new QgsRasterCalcNode( &m ) , 0 );
203+
204+
QgsRasterMatrix result;
205+
result.setNodataValue( -9999 );
206+
QMap<QString, QgsRasterMatrix*> rasterData;
207+
208+
QVERIFY( node.calculate( rasterData, result ) );
209+
210+
QCOMPARE( result.data()[0], -d[0] );
211+
QCOMPARE( result.data()[1], -d[1] );
212+
QCOMPARE( result.data()[2], -d[2] );
213+
QCOMPARE( result.data()[3], -d[3] );
214+
QCOMPARE( result.data()[4], -d[4] );
215+
QCOMPARE( result.data()[5], -9999.0 );
216+
}
217+
218+
void TestQgsRasterCalculator::dualOpNumberMatrix()
219+
{
220+
// test dual op run on number and matrix
221+
double* d = new double[6];
222+
d[0] = 1.0;
223+
d[1] = 2.0;
224+
d[2] = 3.0;
225+
d[3] = 4.0;
226+
d[4] = 5.0;
227+
d[5] = -1.0;
228+
229+
QgsRasterMatrix m( 2, 3, d, -1.0 );
230+
231+
QgsRasterCalcNode node( QgsRasterCalcNode::opPLUS, new QgsRasterCalcNode( 5.0 ), new QgsRasterCalcNode( &m ) );
232+
233+
QgsRasterMatrix result;
234+
result.setNodataValue( -9999 );
235+
QMap<QString, QgsRasterMatrix*> rasterData;
236+
237+
QVERIFY( node.calculate( rasterData, result ) );
238+
239+
QCOMPARE( result.data()[0], 6.0 );
240+
QCOMPARE( result.data()[1], 7.0 );
241+
QCOMPARE( result.data()[2], 8.0 );
242+
QCOMPARE( result.data()[3], 9.0 );
243+
QCOMPARE( result.data()[4], 10.0 );
244+
QCOMPARE( result.data()[5], -9999.0 );
245+
246+
//also check adding no data number
247+
QgsRasterCalcNode nodeNoData( QgsRasterCalcNode::opPLUS, new QgsRasterCalcNode( -9999 ), new QgsRasterCalcNode( &m ) );
248+
QVERIFY( nodeNoData.calculate( rasterData, result ) );
249+
QCOMPARE( result.data()[0], -9999.0 );
250+
}
251+
252+
void TestQgsRasterCalculator::dualOpMatrixNumber()
253+
{
254+
// test dual op run on matrix and number
255+
double* d = new double[6];
256+
d[0] = 1.0;
257+
d[1] = 2.0;
258+
d[2] = 3.0;
259+
d[3] = 4.0;
260+
d[4] = 5.0;
261+
d[5] = -1.0;
262+
263+
QgsRasterMatrix m( 2, 3, d, -1.0 );
264+
265+
QgsRasterCalcNode node( QgsRasterCalcNode::opPLUS, new QgsRasterCalcNode( &m ), new QgsRasterCalcNode( 5.0 ) );
266+
267+
QgsRasterMatrix result;
268+
result.setNodataValue( -9999 );
269+
QMap<QString, QgsRasterMatrix*> rasterData;
270+
271+
QVERIFY( node.calculate( rasterData, result ) );
272+
273+
QCOMPARE( result.data()[0], 6.0 );
274+
QCOMPARE( result.data()[1], 7.0 );
275+
QCOMPARE( result.data()[2], 8.0 );
276+
QCOMPARE( result.data()[3], 9.0 );
277+
QCOMPARE( result.data()[4], 10.0 );
278+
QCOMPARE( result.data()[5], -9999.0 );
279+
280+
//also check adding no data number
281+
QgsRasterCalcNode nodeNoData( QgsRasterCalcNode::opPLUS, new QgsRasterCalcNode( &m ), new QgsRasterCalcNode( -9999 ) );
282+
QVERIFY( nodeNoData.calculate( rasterData, result ) );
283+
QCOMPARE( result.data()[0], -9999.0 );
284+
}
285+
286+
void TestQgsRasterCalculator::dualOpMatrixMatrix()
287+
{
288+
// test dual op run on matrix and matrix
289+
double* d = new double[6];
290+
d[0] = 1.0;
291+
d[1] = 2.0;
292+
d[2] = -2.0;
293+
d[3] = -1.0; //nodata
294+
d[4] = 5.0;
295+
d[5] = -1.0; //nodata
296+
QgsRasterMatrix m1( 2, 3, d, -1.0 );
297+
298+
double* d2 = new double[6];
299+
d2[0] = -1.0;
300+
d2[1] = -2.0; //nodata
301+
d2[2] = 13.0;
302+
d2[3] = -2.0; //nodata
303+
d2[4] = 15.0;
304+
d2[5] = -1.0;
305+
QgsRasterMatrix m2( 2, 3, d2, -2.0 ); //different no data value
306+
307+
QgsRasterCalcNode node( QgsRasterCalcNode::opPLUS, new QgsRasterCalcNode( &m1 ), new QgsRasterCalcNode( &m2 ) );
308+
309+
QgsRasterMatrix result;
310+
result.setNodataValue( -9999 );
311+
QMap<QString, QgsRasterMatrix*> rasterData;
312+
313+
QVERIFY( node.calculate( rasterData, result ) );
314+
315+
QCOMPARE( result.data()[0], 0.0 );
316+
QCOMPARE( result.data()[1], -9999.0 );
317+
QCOMPARE( result.data()[2], 11.0 );
318+
QCOMPARE( result.data()[3], -9999.0 );
319+
QCOMPARE( result.data()[4], 20.0 );
320+
QCOMPARE( result.data()[5], -9999.0 );
321+
}
322+
323+
void TestQgsRasterCalculator::rasterRefOp()
324+
{
325+
// test single op run on raster ref
326+
QgsRasterCalcNode node( QgsRasterCalcNode::opSIGN, new QgsRasterCalcNode( "raster" ), 0 );
327+
328+
QgsRasterMatrix result;
329+
result.setNodataValue( -9999 );
330+
QMap<QString, QgsRasterMatrix*> rasterData;
331+
332+
//first test invalid raster ref
333+
QVERIFY( !node.calculate( rasterData, result ) );
334+
335+
//now create raster ref
336+
double* d = new double[6];
337+
d[0] = 1.0;
338+
d[1] = 2.0;
339+
d[2] = 3.0;
340+
d[3] = 4.0;
341+
d[4] = 5.0;
342+
d[5] = -1.0;
343+
QgsRasterMatrix m( 2, 3, d, -1.0 );
344+
rasterData.insert( "raster", &m );
345+
346+
QVERIFY( node.calculate( rasterData, result ) );
347+
QCOMPARE( result.data()[0], -d[0] );
348+
QCOMPARE( result.data()[1], -d[1] );
349+
QCOMPARE( result.data()[2], -d[2] );
350+
QCOMPARE( result.data()[3], -d[3] );
351+
QCOMPARE( result.data()[4], -d[4] );
352+
QCOMPARE( result.data()[5], -9999.0 );
353+
}
354+
355+
void TestQgsRasterCalculator::dualOpRasterRaster()
356+
{
357+
// test dual op run on matrix and matrix
358+
double* d = new double[6];
359+
d[0] = 1.0;
360+
d[1] = 2.0;
361+
d[2] = -2.0;
362+
d[3] = -1.0; //nodata
363+
d[4] = 5.0;
364+
d[5] = -1.0; //nodata
365+
QgsRasterMatrix m1( 2, 3, d, -1.0 );
366+
QMap<QString, QgsRasterMatrix*> rasterData;
367+
rasterData.insert( "raster1", &m1 );
368+
369+
double* d2 = new double[6];
370+
d2[0] = -1.0;
371+
d2[1] = -2.0; //nodata
372+
d2[2] = 13.0;
373+
d2[3] = -2.0; //nodata
374+
d2[4] = 15.0;
375+
d2[5] = -1.0;
376+
QgsRasterMatrix m2( 2, 3, d2, -2.0 ); //different no data value
377+
rasterData.insert( "raster2", &m2 );
378+
379+
QgsRasterCalcNode node( QgsRasterCalcNode::opPLUS, new QgsRasterCalcNode( "raster1" ), new QgsRasterCalcNode( "raster2" ) );
380+
381+
QgsRasterMatrix result;
382+
result.setNodataValue( -9999 );
383+
384+
QVERIFY( node.calculate( rasterData, result ) );
385+
QCOMPARE( result.data()[0], 0.0 );
386+
QCOMPARE( result.data()[1], -9999.0 );
387+
QCOMPARE( result.data()[2], 11.0 );
388+
QCOMPARE( result.data()[3], -9999.0 );
389+
QCOMPARE( result.data()[4], 20.0 );
390+
QCOMPARE( result.data()[5], -9999.0 );
391+
}
392+
181393
QTEST_MAIN( TestQgsRasterCalculator )
182394
#include "testqgsrastercalculator.moc"

0 commit comments

Comments
 (0)
Please sign in to comment.