QGIS API Documentation 4.1.0-Master (d6fb7a379fb)
Loading...
Searching...
No Matches
qgsrastercalcnode.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsrastercalcnode.cpp
3 ---------------------
4 begin : October 2010
5 copyright : (C) 2010 by Marco Hugentobler
6 email : marco dot hugentobler at sourcepole dot ch
7 ***************************************************************************
8 * *
9 * This program is free software; you can redistribute it and/or modify *
10 * it under the terms of the GNU General Public License as published by *
11 * the Free Software Foundation; either version 2 of the License, or *
12 * (at your option) any later version. *
13 * *
14 ***************************************************************************/
15#include "qgsrastercalcnode.h"
16
17#include "qgsrasterblock.h"
18#include "qgsrastermatrix.h"
19
20#include <QString>
21
22using namespace Qt::StringLiterals;
23
25 : mNumber( number )
26{}
27
29 : mType( tMatrix )
30 , mMatrix( matrix )
31{}
32
34 : mType( tOperator )
35 , mLeft( left )
36 , mRight( right )
37 , mOperator( op )
38{}
39
40QgsRasterCalcNode::QgsRasterCalcNode( QString functionName, QVector<QgsRasterCalcNode *> functionArgs )
41 : mType( tFunction )
42 , mFunctionName( functionName )
43 , mFunctionArgs( functionArgs )
44{}
45
46QgsRasterCalcNode::QgsRasterCalcNode( const QString &rasterName )
47 : mType( tRasterRef )
48 , mRasterName( rasterName )
49{
50 if ( mRasterName.startsWith( '"' ) && mRasterName.endsWith( '"' ) )
51 mRasterName = mRasterName.mid( 1, mRasterName.size() - 2 );
52}
53
55{
56 for ( int i = 0; i < mFunctionArgs.size(); ++i )
57 {
58 if ( mFunctionArgs.at( i ) )
59 delete mFunctionArgs.at( i );
60 }
61}
62
63bool QgsRasterCalcNode::calculate( QMap<QString, QgsRasterBlock *> &rasterData, QgsRasterMatrix &result, int row ) const
64{
65 //if type is raster ref: return a copy of the corresponding matrix
66
67 //if type is operator, call the proper matrix operations
68 if ( mType == tRasterRef )
69 {
70 const QMap<QString, QgsRasterBlock *>::iterator it = rasterData.find( mRasterName );
71 if ( it == rasterData.end() )
72 {
73 QgsDebugError( u"Error: could not find raster data for \"%1\""_s.arg( mRasterName ) );
74 return false;
75 }
76
77 const int nRows = ( row >= 0 ? 1 : ( *it )->height() );
78 const int startRow = ( row >= 0 ? row : 0 );
79 const int endRow = startRow + nRows;
80 const int nCols = ( *it )->width();
81 const int nEntries = nCols * nRows;
82 double *data = new double[nEntries];
83
84 //convert input raster values to double, also convert input no data to result no data
85
86 int outRow = 0;
87 bool isNoData = false;
88 for ( int dataRow = startRow; dataRow < endRow; ++dataRow, ++outRow )
89 {
90 for ( int dataCol = 0; dataCol < nCols; ++dataCol )
91 {
92 const double value = ( *it )->valueAndNoData( dataRow, dataCol, isNoData );
93 data[dataCol + nCols * outRow] = isNoData ? result.nodataValue() : value;
94 }
95 }
96 result.setData( nCols, nRows, data, result.nodataValue() );
97 return true;
98 }
99 else if ( mType == tOperator )
100 {
101 QgsRasterMatrix leftMatrix( result.nColumns(), result.nRows(), nullptr, result.nodataValue() );
102 QgsRasterMatrix rightMatrix( result.nColumns(), result.nRows(), nullptr, result.nodataValue() );
103
104 if ( !mLeft || !mLeft->calculate( rasterData, leftMatrix, row ) )
105 {
106 return false;
107 }
108 if ( mRight && !mRight->calculate( rasterData, rightMatrix, row ) )
109 {
110 return false;
111 }
112
113 switch ( mOperator )
114 {
115 case opPLUS:
116 leftMatrix.add( rightMatrix );
117 break;
118 case opMINUS:
119 leftMatrix.subtract( rightMatrix );
120 break;
121 case opMUL:
122 leftMatrix.multiply( rightMatrix );
123 break;
124 case opDIV:
125 leftMatrix.divide( rightMatrix );
126 break;
127 case opPOW:
128 leftMatrix.power( rightMatrix );
129 break;
130 case opEQ:
131 leftMatrix.equal( rightMatrix );
132 break;
133 case opNE:
134 leftMatrix.notEqual( rightMatrix );
135 break;
136 case opGT:
137 leftMatrix.greaterThan( rightMatrix );
138 break;
139 case opLT:
140 leftMatrix.lesserThan( rightMatrix );
141 break;
142 case opGE:
143 leftMatrix.greaterEqual( rightMatrix );
144 break;
145 case opLE:
146 leftMatrix.lesserEqual( rightMatrix );
147 break;
148 case opAND:
149 leftMatrix.logicalAnd( rightMatrix );
150 break;
151 case opOR:
152 leftMatrix.logicalOr( rightMatrix );
153 break;
154 case opMIN:
155 leftMatrix.min( rightMatrix );
156 break;
157 case opMAX:
158 leftMatrix.max( rightMatrix );
159 break;
160 case opSQRT:
161 leftMatrix.squareRoot();
162 break;
163 case opSIN:
164 leftMatrix.sinus();
165 break;
166 case opCOS:
167 leftMatrix.cosinus();
168 break;
169 case opTAN:
170 leftMatrix.tangens();
171 break;
172 case opASIN:
173 leftMatrix.asinus();
174 break;
175 case opACOS:
176 leftMatrix.acosinus();
177 break;
178 case opATAN:
179 leftMatrix.atangens();
180 break;
181 case opSIGN:
182 leftMatrix.changeSign();
183 break;
184 case opLOG:
185 leftMatrix.log();
186 break;
187 case opLOG10:
188 leftMatrix.log10();
189 break;
190 case opABS:
191 leftMatrix.absoluteValue();
192 break;
193 default:
194 return false;
195 }
196 const int newNColumns = leftMatrix.nColumns();
197 const int newNRows = leftMatrix.nRows();
198 result.setData( newNColumns, newNRows, leftMatrix.takeData(), leftMatrix.nodataValue() );
199 return true;
200 }
201 else if ( mType == tNumber )
202 {
203 const size_t nEntries = static_cast<size_t>( result.nColumns() ) * static_cast<size_t>( result.nRows() );
204 double *data = new double[nEntries];
205 std::fill( data, data + nEntries, mNumber );
206 result.setData( result.nColumns(), 1, data, result.nodataValue() );
207
208 return true;
209 }
210 else if ( mType == tMatrix )
211 {
212 const int nEntries = mMatrix->nColumns() * mMatrix->nRows();
213 double *data = new double[nEntries];
214 for ( int i = 0; i < nEntries; ++i )
215 {
216 data[i] = mMatrix->data()[i] == mMatrix->nodataValue() ? result.nodataValue() : mMatrix->data()[i];
217 }
218 result.setData( mMatrix->nColumns(), mMatrix->nRows(), data, result.nodataValue() );
219 return true;
220 }
221 else if ( mType == tFunction )
222 {
223 std::vector<std::unique_ptr<QgsRasterMatrix>> matrixContainer;
224 for ( int i = 0; i < mFunctionArgs.size(); ++i )
225 {
226 auto singleMatrix = std::make_unique<QgsRasterMatrix>( result.nColumns(), result.nRows(), nullptr, result.nodataValue() );
227 if ( !mFunctionArgs.at( i ) || !mFunctionArgs.at( i )->calculate( rasterData, *singleMatrix, row ) )
228 {
229 return false;
230 }
231 matrixContainer.emplace_back( std::move( singleMatrix ) );
232 }
233 evaluateFunction( matrixContainer, result );
234 return true;
235 }
236 return false;
237}
238
239QString QgsRasterCalcNode::toString( bool cStyle ) const
240{
241 QString result;
242 QString left;
243 QString right;
244 if ( mLeft )
245 left = mLeft->toString( cStyle );
246 if ( mRight )
247 right = mRight->toString( cStyle );
248
249 switch ( mType )
250 {
251 case tOperator:
252 switch ( mOperator )
253 {
254 case opPLUS:
255 result = u"( %1 + %2 )"_s.arg( left, right );
256 break;
257 case opMINUS:
258 result = u"( %1 - %2 )"_s.arg( left, right );
259 break;
260 case opSIGN:
261 result = u"-%1"_s.arg( left );
262 break;
263 case opMUL:
264 result = u"%1 * %2"_s.arg( left, right );
265 break;
266 case opDIV:
267 result = u"%1 / %2"_s.arg( left, right );
268 break;
269 case opPOW:
270 if ( cStyle )
271 result = u"pow( %1, %2 )"_s.arg( left, right );
272 else
273 result = u"%1^%2"_s.arg( left, right );
274 break;
275 case opEQ:
276 if ( cStyle )
277 result = u"( float ) ( %1 == %2 )"_s.arg( left, right );
278 else
279 result = u"%1 = %2"_s.arg( left, right );
280 break;
281 case opNE:
282 if ( cStyle )
283 result = u"( float ) ( %1 != %2 )"_s.arg( left, right );
284 else
285 result = u"%1 != %2"_s.arg( left, right );
286 break;
287 case opGT:
288 if ( cStyle )
289 result = u"( float ) ( %1 > %2 )"_s.arg( left, right );
290 else
291 result = u"%1 > %2"_s.arg( left, right );
292 break;
293 case opLT:
294 if ( cStyle )
295 result = u"( float ) ( %1 < %2 )"_s.arg( left, right );
296 else
297 result = u"%1 < %2"_s.arg( left, right );
298 break;
299 case opGE:
300 if ( cStyle )
301 result = u"( float ) ( %1 >= %2 )"_s.arg( left, right );
302 else
303 result = u"%1 >= %2"_s.arg( left, right );
304 break;
305 case opLE:
306 if ( cStyle )
307 result = u"( float ) ( %1 <= %2 )"_s.arg( left, right );
308 else
309 result = u"%1 <= %2"_s.arg( left, right );
310 break;
311 case opAND:
312 if ( cStyle )
313 result = u"( float ) ( %1 && %2 )"_s.arg( left, right );
314 else
315 result = u"%1 AND %2"_s.arg( left, right );
316 break;
317 case opOR:
318 if ( cStyle )
319 result = u"( float ) ( %1 || %2 )"_s.arg( left, right );
320 else
321 result = u"%1 OR %2"_s.arg( left, right );
322 break;
323 case opSQRT:
324 result = u"sqrt( %1 )"_s.arg( left );
325 break;
326 case opSIN:
327 result = u"sin( %1 )"_s.arg( left );
328 break;
329 case opCOS:
330 result = u"cos( %1 )"_s.arg( left );
331 break;
332 case opTAN:
333 result = u"tan( %1 )"_s.arg( left );
334 break;
335 case opASIN:
336 result = u"asin( %1 )"_s.arg( left );
337 break;
338 case opACOS:
339 result = u"acos( %1 )"_s.arg( left );
340 break;
341 case opATAN:
342 result = u"atan( %1 )"_s.arg( left );
343 break;
344 case opLOG:
345 result = u"log( %1 )"_s.arg( left );
346 break;
347 case opLOG10:
348 result = u"log10( %1 )"_s.arg( left );
349 break;
350 case opABS:
351 if ( cStyle )
352 result = u"fabs( %1 )"_s.arg( left );
353 else
354 // Call the floating point version
355 result = u"abs( %1 )"_s.arg( left );
356 break;
357 case opMIN:
358 if ( cStyle )
359 result = u"min( ( float ) ( %1 ), ( float ) ( %2 ) )"_s.arg( left, right );
360 else
361 result = u"min( %1, %2 )"_s.arg( left, right );
362 break;
363 case opMAX:
364 if ( cStyle )
365 result = u"max( ( float ) ( %1 ), ( float ) ( %2 ) )"_s.arg( left, right );
366 else
367 result = u"max( %1, %2 )"_s.arg( left, right );
368 break;
369 case opNONE:
370 break;
371 }
372 break;
373 case tRasterRef:
374 if ( cStyle )
375 result = u"( float ) \"%1\""_s.arg( mRasterName );
376 else
377 result = u"\"%1\""_s.arg( mRasterName );
378 break;
379 case tNumber:
380 result = QString::number( mNumber );
381 if ( cStyle )
382 {
383 result = u"( float ) %1"_s.arg( result );
384 }
385 break;
386 case tMatrix:
387 break;
388 case tFunction:
389 if ( mFunctionName == "if" )
390 {
391 const QString argOne = mFunctionArgs.at( 0 )->toString( cStyle );
392 const QString argTwo = mFunctionArgs.at( 1 )->toString( cStyle );
393 const QString argThree = mFunctionArgs.at( 2 )->toString( cStyle );
394 if ( cStyle )
395 result = u" ( bool ) ( %1 ) ? ( %2 ) : ( %3 ) "_s.arg( argOne, argTwo, argThree );
396 else
397 result = u"if( %1 , %2 , %3 )"_s.arg( argOne, argTwo, argThree );
398 }
399 break;
400 }
401 return result;
402}
403
404QList<const QgsRasterCalcNode *> QgsRasterCalcNode::findNodes( const QgsRasterCalcNode::Type type ) const
405{
406 QList<const QgsRasterCalcNode *> nodeList;
407 if ( mType == type )
408 nodeList.push_back( this );
409 if ( mLeft )
410 nodeList.append( mLeft->findNodes( type ) );
411 if ( mRight )
412 nodeList.append( mRight->findNodes( type ) );
413
414 for ( QgsRasterCalcNode *node : mFunctionArgs )
415 nodeList.append( node->findNodes( type ) );
416
417 return nodeList;
418}
419
420QgsRasterCalcNode *QgsRasterCalcNode::parseRasterCalcString( const QString &str, QString &parserErrorMsg )
421{
422 extern QgsRasterCalcNode *localParseRasterCalcString( const QString &str, QString &parserErrorMsg );
423 return localParseRasterCalcString( str, parserErrorMsg );
424}
425
427{
428 QStringList referencedRasters;
429
430 const QStringList rasterRef = this->cleanRasterReferences();
431 for ( const auto &i : std::as_const( rasterRef ) )
432 {
433 if ( referencedRasters.contains( i.mid( 0, i.lastIndexOf( "@" ) ) ) )
434 continue;
435 referencedRasters << i.mid( 0, i.lastIndexOf( "@" ) );
436 }
437
438 return referencedRasters;
439}
440
442{
443 QStringList rasterReferences;
444 const QList<const QgsRasterCalcNode *> rasterRefNodes = this->findNodes( QgsRasterCalcNode::Type::tRasterRef );
445
446 for ( const QgsRasterCalcNode *r : rasterRefNodes )
447 {
448 QString layerRef( r->toString() );
449 if ( layerRef.at( 0 ) == "\""_L1 && layerRef.at( layerRef.size() - 1 ) == "\""_L1 )
450 {
451 layerRef.remove( 0, 1 );
452 layerRef.chop( 1 );
453 }
454 layerRef.remove( QChar( '\\' ), Qt::CaseInsensitive );
455 rasterReferences << layerRef;
456 }
457
458 return rasterReferences;
459}
460
461QgsRasterMatrix QgsRasterCalcNode::evaluateFunction( const std::vector<std::unique_ptr<QgsRasterMatrix>> &matrixVector, QgsRasterMatrix &result ) const
462{
463 if ( mFunctionName == "if" )
464 {
465 //scalar condition
466 if ( matrixVector.at( 0 )->isNumber() )
467 {
468 result = ( matrixVector.at( 0 )->data() ? *matrixVector.at( 1 ) : *matrixVector.at( 2 ) );
469 return result;
470 }
471 int nCols = matrixVector.at( 0 )->nColumns();
472 int nRows = matrixVector.at( 0 )->nRows();
473 int nEntries = nCols * nRows;
474 auto dataResult = std::make_unique<double[]>( nEntries );
475 double *dataResultRawPtr = dataResult.get();
476
477 double *condition = matrixVector.at( 0 )->data();
478 double *firstOption = matrixVector.at( 1 )->data();
479 double *secondOption = matrixVector.at( 2 )->data();
480
481 bool isFirstOptionNumber = matrixVector.at( 1 )->isNumber();
482 bool isSecondCOptionNumber = matrixVector.at( 2 )->isNumber();
483 double noDataValueCondition = matrixVector.at( 0 )->nodataValue();
484
485 for ( int i = 0; i < nEntries; ++i )
486 {
487 if ( condition[i] == noDataValueCondition )
488 {
489 dataResultRawPtr[i] = result.nodataValue();
490 continue;
491 }
492 else if ( condition[i] != 0 )
493 {
494 dataResultRawPtr[i] = isFirstOptionNumber ? firstOption[0] : firstOption[i];
495 continue;
496 }
497 dataResultRawPtr[i] = isSecondCOptionNumber ? secondOption[0] : secondOption[i];
498 }
499
500 result.setData( nCols, nRows, dataResult.release(), result.nodataValue() );
501 }
502 return result;
503}
QString toString(bool cStyle=false) const
Returns a string representation of the expression.
Operator
possible operators
QgsRasterCalcNode()=default
QList< const QgsRasterCalcNode * > findNodes(const QgsRasterCalcNode::Type type) const
Returns a list of nodes of a specific type.
static QgsRasterCalcNode * parseRasterCalcString(const QString &str, QString &parserErrorMsg)
bool calculate(QMap< QString, QgsRasterBlock * > &rasterData, QgsRasterMatrix &result, int row=-1) const
Calculates result of raster calculation (might be real matrix or single number).
QStringList referencedLayerNames() const
Returns a list of raster layer names that are referenced in the formula without the quotation marks.
QStringList cleanRasterReferences() const
Returns a list of raster layer references that are addressed in the formula, without quotation marks.
Type
defines possible types of node
Represents a matrix in a raster calculator operation.
bool equal(const QgsRasterMatrix &other)
bool max(const QgsRasterMatrix &other)
Calculates the maximum value between two matrices.
bool greaterEqual(const QgsRasterMatrix &other)
double * takeData()
Returns data and ownership.
bool notEqual(const QgsRasterMatrix &other)
bool greaterThan(const QgsRasterMatrix &other)
bool logicalOr(const QgsRasterMatrix &other)
bool lesserEqual(const QgsRasterMatrix &other)
bool absoluteValue()
Calculates the absolute value.
bool logicalAnd(const QgsRasterMatrix &other)
bool lesserThan(const QgsRasterMatrix &other)
bool add(const QgsRasterMatrix &other)
Adds another matrix to this one.
bool multiply(const QgsRasterMatrix &other)
int nRows() const
bool power(const QgsRasterMatrix &other)
double nodataValue() const
bool min(const QgsRasterMatrix &other)
Calculates the minimum value between two matrices.
bool divide(const QgsRasterMatrix &other)
int nColumns() const
void setData(int cols, int rows, double *data, double nodataValue)
bool subtract(const QgsRasterMatrix &other)
Subtracts another matrix from this one.
#define QgsDebugError(str)
Definition qgslogger.h:59