Statistics
| Branch: | Tag: | Revision:

qgis / src / analysis / raster / qgsrastercalcnode.cpp @ master

History | View | Annotate | Download (15.8 KB)

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
#include "qgsrasterblock.h"
17
#include "qgsrastermatrix.h"
18

    
19
QgsRasterCalcNode::QgsRasterCalcNode( double number )
20
  : mNumber( number )
21
{
22
}
23

    
24
QgsRasterCalcNode::QgsRasterCalcNode( QgsRasterMatrix *matrix )
25
  : mType( tMatrix )
26
  , mMatrix( matrix )
27
{
28

    
29
}
30

    
31
QgsRasterCalcNode::QgsRasterCalcNode( Operator op, QgsRasterCalcNode *left, QgsRasterCalcNode *right )
32
  : mType( tOperator )
33
  , mLeft( left )
34
  , mRight( right )
35
  , mOperator( op )
36
{
37
}
38

    
39
QgsRasterCalcNode::QgsRasterCalcNode( QString functionName, QVector <QgsRasterCalcNode *> functionArgs )
40
  : mType( tFunction )
41
  , mFunctionName( functionName )
42
  , mFunctionArgs( functionArgs )
43
{
44
}
45

    
46
QgsRasterCalcNode::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

    
54
QgsRasterCalcNode::~QgsRasterCalcNode()
55
{
56
  delete mLeft;
57
  delete mRight;
58
}
59

    
60
bool QgsRasterCalcNode::calculate( QMap<QString, QgsRasterBlock * > &rasterData, QgsRasterMatrix &result, int row ) const
61
{
62
  //if type is raster ref: return a copy of the corresponding matrix
63

    
64
  //if type is operator, call the proper matrix operations
65
  if ( mType == tRasterRef )
66
  {
67
    const QMap<QString, QgsRasterBlock *>::iterator it = rasterData.find( mRasterName );
68
    if ( it == rasterData.end() )
69
    {
70
      QgsDebugMsg( QStringLiteral( "Error: could not find raster data for \"%1\"" ).arg( mRasterName ) );
71
      return false;
72
    }
73

    
74
    const int nRows = ( row >= 0 ? 1 : ( *it )->height() );
75
    const int startRow = ( row >= 0 ? row : 0 );
76
    const int endRow = startRow + nRows;
77
    const int nCols = ( *it )->width();
78
    const int nEntries = nCols * nRows;
79
    double *data = new double[nEntries];
80

    
81
    //convert input raster values to double, also convert input no data to result no data
82

    
83
    int outRow = 0;
84
    bool isNoData = false;
85
    for ( int dataRow = startRow; dataRow < endRow ; ++dataRow, ++outRow )
86
    {
87
      for ( int dataCol = 0; dataCol < nCols; ++dataCol )
88
      {
89
        const double value = ( *it )->valueAndNoData( dataRow, dataCol, isNoData );
90
        data[ dataCol + nCols * outRow] = isNoData ? result.nodataValue() : value;
91
      }
92
    }
93
    result.setData( nCols, nRows, data, result.nodataValue() );
94
    return true;
95
  }
96
  else if ( mType == tOperator )
97
  {
98
    QgsRasterMatrix leftMatrix( result.nColumns(), result.nRows(), nullptr, result.nodataValue() );
99
    QgsRasterMatrix rightMatrix( result.nColumns(), result.nRows(), nullptr, result.nodataValue() );
100

    
101
    if ( !mLeft || !mLeft->calculate( rasterData, leftMatrix, row ) )
102
    {
103
      return false;
104
    }
105
    if ( mRight && !mRight->calculate( rasterData, rightMatrix, row ) )
106
    {
107
      return false;
108
    }
109

    
110
    switch ( mOperator )
111
    {
112
      case opPLUS:
113
        leftMatrix.add( rightMatrix );
114
        break;
115
      case opMINUS:
116
        leftMatrix.subtract( rightMatrix );
117
        break;
118
      case opMUL:
119
        leftMatrix.multiply( rightMatrix );
120
        break;
121
      case opDIV:
122
        leftMatrix.divide( rightMatrix );
123
        break;
124
      case opPOW:
125
        leftMatrix.power( rightMatrix );
126
        break;
127
      case opEQ:
128
        leftMatrix.equal( rightMatrix );
129
        break;
130
      case opNE:
131
        leftMatrix.notEqual( rightMatrix );
132
        break;
133
      case opGT:
134
        leftMatrix.greaterThan( rightMatrix );
135
        break;
136
      case opLT:
137
        leftMatrix.lesserThan( rightMatrix );
138
        break;
139
      case opGE:
140
        leftMatrix.greaterEqual( rightMatrix );
141
        break;
142
      case opLE:
143
        leftMatrix.lesserEqual( rightMatrix );
144
        break;
145
      case opAND:
146
        leftMatrix.logicalAnd( rightMatrix );
147
        break;
148
      case opOR:
149
        leftMatrix.logicalOr( rightMatrix );
150
        break;
151
      case opMIN:
152
        leftMatrix.min( rightMatrix );
153
        break;
154
      case opMAX:
155
        leftMatrix.max( rightMatrix );
156
        break;
157
      case opSQRT:
158
        leftMatrix.squareRoot();
159
        break;
160
      case opSIN:
161
        leftMatrix.sinus();
162
        break;
163
      case opCOS:
164
        leftMatrix.cosinus();
165
        break;
166
      case opTAN:
167
        leftMatrix.tangens();
168
        break;
169
      case opASIN:
170
        leftMatrix.asinus();
171
        break;
172
      case opACOS:
173
        leftMatrix.acosinus();
174
        break;
175
      case opATAN:
176
        leftMatrix.atangens();
177
        break;
178
      case opSIGN:
179
        leftMatrix.changeSign();
180
        break;
181
      case opLOG:
182
        leftMatrix.log();
183
        break;
184
      case opLOG10:
185
        leftMatrix.log10();
186
        break;
187
      case opABS:
188
        leftMatrix.absoluteValue();
189
        break;
190
      default:
191
        return false;
192
    }
193
    const int newNColumns = leftMatrix.nColumns();
194
    const int newNRows = leftMatrix.nRows();
195
    result.setData( newNColumns, newNRows, leftMatrix.takeData(), leftMatrix.nodataValue() );
196
    return true;
197
  }
198
  else if ( mType == tNumber )
199
  {
200
    const size_t nEntries = static_cast<size_t>( result.nColumns() * result.nRows() );
201
    double *data = new double[ nEntries ];
202
    std::fill( data, data + nEntries, mNumber );
203
    result.setData( result.nColumns(), 1, data, result.nodataValue() );
204

    
205
    return true;
206
  }
207
  else if ( mType == tMatrix )
208
  {
209
    const int nEntries = mMatrix->nColumns() * mMatrix->nRows();
210
    double *data = new double[nEntries];
211
    for ( int i = 0; i < nEntries; ++i )
212
    {
213
      data[i] = mMatrix->data()[i] == mMatrix->nodataValue() ? result.nodataValue() : mMatrix->data()[i];
214
    }
215
    result.setData( mMatrix->nColumns(), mMatrix->nRows(), data, result.nodataValue() );
216
    return true;
217
  }
218
  else if ( mType == tFunction )
219
  {
220
    QVector <QgsRasterMatrix *> matrixContainer;
221
    for ( int i = 0; i < mFunctionArgs.size(); ++i )
222
    {
223
      std::unique_ptr< QgsRasterMatrix > singleMatrix( new QgsRasterMatrix( result.nColumns(), result.nRows(), nullptr, result.nodataValue() ) );
224
      if ( !mFunctionArgs.at( i ) || !mFunctionArgs.at( i )->calculate( rasterData, *singleMatrix, row ) )
225
      {
226
        return false;
227
      }
228
      matrixContainer.append( singleMatrix.release() );
229
    }
230
    evaluateFunction( matrixContainer, result );
231
    return true;
232
  }
233
  return false;
234
}
235

    
236
QString QgsRasterCalcNode::toString( bool cStyle ) const
237
{
238
  QString result;
239
  QString left;
240
  QString right;
241
  if ( mLeft )
242
    left = mLeft->toString( cStyle );
243
  if ( mRight )
244
    right = mRight->toString( cStyle );
245

    
246
  switch ( mType )
247
  {
248
    case tOperator:
249
      switch ( mOperator )
250
      {
251
        case opPLUS:
252
          result = QStringLiteral( "( %1 + %2 )" ).arg( left ).arg( right );
253
          break;
254
        case opMINUS:
255
          result = QStringLiteral( "( %1 - %2 )" ).arg( left ).arg( right );
256
          break;
257
        case opSIGN:
258
          result = QStringLiteral( "-%1" ).arg( left );
259
          break;
260
        case opMUL:
261
          result = QStringLiteral( "%1 * %2" ).arg( left ).arg( right );
262
          break;
263
        case opDIV:
264
          result = QStringLiteral( "%1 / %2" ).arg( left ).arg( right );
265
          break;
266
        case opPOW:
267
          if ( cStyle )
268
            result = QStringLiteral( "pow( %1, %2 )" ).arg( left ).arg( right );
269
          else
270
            result = QStringLiteral( "%1^%2" ).arg( left ).arg( right );
271
          break;
272
        case opEQ:
273
          if ( cStyle )
274
            result = QStringLiteral( "( float ) ( %1 == %2 )" ).arg( left ).arg( right );
275
          else
276
            result = QStringLiteral( "%1 = %2" ).arg( left ).arg( right );
277
          break;
278
        case opNE:
279
          if ( cStyle )
280
            result = QStringLiteral( "( float ) ( %1 != %2 )" ).arg( left ).arg( right );
281
          else
282
            result = QStringLiteral( "%1 != %2" ).arg( left ).arg( right );
283
          break;
284
        case opGT:
285
          if ( cStyle )
286
            result = QStringLiteral( "( float ) ( %1 > %2 )" ).arg( left ).arg( right );
287
          else
288
            result = QStringLiteral( "%1 > %2" ).arg( left ).arg( right );
289
          break;
290
        case opLT:
291
          if ( cStyle )
292
            result = QStringLiteral( "( float ) ( %1 < %2 )" ).arg( left ).arg( right );
293
          else
294
            result = QStringLiteral( "%1 < %2" ).arg( left ).arg( right );
295
          break;
296
        case opGE:
297
          if ( cStyle )
298
            result = QStringLiteral( "( float ) ( %1 >= %2 )" ).arg( left ).arg( right );
299
          else
300
            result = QStringLiteral( "%1 >= %2" ).arg( left ).arg( right );
301
          break;
302
        case opLE:
303
          if ( cStyle )
304
            result = QStringLiteral( "( float ) ( %1 <= %2 )" ).arg( left ).arg( right );
305
          else
306
            result = QStringLiteral( "%1 <= %2" ).arg( left ).arg( right );
307
          break;
308
        case opAND:
309
          if ( cStyle )
310
            result = QStringLiteral( "( float ) ( %1 && %2 )" ).arg( left ).arg( right );
311
          else
312
            result = QStringLiteral( "%1 AND %2" ).arg( left ).arg( right );
313
          break;
314
        case opOR:
315
          if ( cStyle )
316
            result = QStringLiteral( "( float ) ( %1 || %2 )" ).arg( left ).arg( right );
317
          else
318
            result = QStringLiteral( "%1 OR %2" ).arg( left ).arg( right );
319
          break;
320
        case opSQRT:
321
          result = QStringLiteral( "sqrt( %1 )" ).arg( left );
322
          break;
323
        case opSIN:
324
          result = QStringLiteral( "sin( %1 )" ).arg( left );
325
          break;
326
        case opCOS:
327
          result = QStringLiteral( "cos( %1 )" ).arg( left );
328
          break;
329
        case opTAN:
330
          result = QStringLiteral( "tan( %1 )" ).arg( left );
331
          break;
332
        case opASIN:
333
          result = QStringLiteral( "asin( %1 )" ).arg( left );
334
          break;
335
        case opACOS:
336
          result = QStringLiteral( "acos( %1 )" ).arg( left );
337
          break;
338
        case opATAN:
339
          result = QStringLiteral( "atan( %1 )" ).arg( left );
340
          break;
341
        case opLOG:
342
          result = QStringLiteral( "log( %1 )" ).arg( left );
343
          break;
344
        case opLOG10:
345
          result = QStringLiteral( "log10( %1 )" ).arg( left );
346
          break;
347
        case opABS:
348
          if ( cStyle )
349
            result = QStringLiteral( "fabs( %1 )" ).arg( left );
350
          else
351
            // Call the floating point version
352
            result = QStringLiteral( "abs( %1 )" ).arg( left );
353
          break;
354
        case opMIN:
355
          if ( cStyle )
356
            result = QStringLiteral( "min( ( float ) ( %1 ), ( float ) ( %2 ) )" ).arg( left ).arg( right );
357
          else
358
            result = QStringLiteral( "min( %1, %2 )" ).arg( left ).arg( right );
359
          break;
360
        case opMAX:
361
          if ( cStyle )
362
            result = QStringLiteral( "max( ( float ) ( %1 ), ( float ) ( %2 ) )" ).arg( left ).arg( right );
363
          else
364
            result = QStringLiteral( "max( %1, %2 )" ).arg( left ).arg( right );
365
          break;
366
        case opNONE:
367
          break;
368
      }
369
      break;
370
    case tRasterRef:
371
      if ( cStyle )
372
        result = QStringLiteral( "( float ) \"%1\"" ).arg( mRasterName );
373
      else
374
        result = QStringLiteral( "\"%1\"" ).arg( mRasterName );
375
      break;
376
    case tNumber:
377
      result = QString::number( mNumber );
378
      if ( cStyle )
379
      {
380
        result = QStringLiteral( "( float ) %1" ).arg( result );
381
      }
382
      break;
383
    case tMatrix:
384
      break;
385
    case tFunction:
386
      if ( mFunctionName == "if" )
387
      {
388
        const QString argOne = mFunctionArgs.at( 0 )->toString( cStyle );
389
        const QString argTwo = mFunctionArgs.at( 1 )->toString( cStyle );
390
        const QString argThree = mFunctionArgs.at( 2 )->toString( cStyle );
391
        if ( cStyle )
392
          result =  QStringLiteral( " ( %1 ) ? ( %2 ) : ( %3 ) " ).arg( argOne, argTwo, argThree );
393
        else
394
          result = QStringLiteral( "if( %1 , %2 , %3 )" ).arg( argOne, argTwo, argThree );
395
      }
396
      break;
397
  }
398
  return result;
399
}
400

    
401
QList<const QgsRasterCalcNode *> QgsRasterCalcNode::findNodes( const QgsRasterCalcNode::Type type ) const
402
{
403
  QList<const QgsRasterCalcNode *> nodeList;
404
  if ( mType == type )
405
    nodeList.push_back( this );
406
  if ( mLeft )
407
    nodeList.append( mLeft->findNodes( type ) );
408
  if ( mRight )
409
    nodeList.append( mRight->findNodes( type ) );
410

    
411
  for ( QgsRasterCalcNode *node : mFunctionArgs )
412
    nodeList.append( node->findNodes( type ) );
413

    
414
  return nodeList;
415
}
416

    
417
QgsRasterCalcNode *QgsRasterCalcNode::parseRasterCalcString( const QString &str, QString &parserErrorMsg )
418
{
419
  extern QgsRasterCalcNode *localParseRasterCalcString( const QString & str, QString & parserErrorMsg );
420
  return localParseRasterCalcString( str, parserErrorMsg );
421
}
422

    
423
QStringList QgsRasterCalcNode::referencedLayerNames()
424
{
425
  QStringList referencedRasters;
426

    
427
  QStringList rasterRef = this->cleanRasterReferences();
428
  for ( const auto &i : rasterRef )
429
  {
430
    if ( referencedRasters.contains( i.mid( 0, i.lastIndexOf( "@" ) ) ) ) continue;
431
    referencedRasters << i.mid( 0, i.lastIndexOf( "@" ) );
432
  }
433

    
434
  return referencedRasters;
435
}
436

    
437
QStringList QgsRasterCalcNode::cleanRasterReferences()
438
{
439
  QStringList rasterReferences;
440
  const QList<const QgsRasterCalcNode *> rasterRefNodes =  this->findNodes( QgsRasterCalcNode::Type::tRasterRef );
441

    
442
  for ( const QgsRasterCalcNode *r : rasterRefNodes )
443
  {
444

    
445
    QString layerRef( r->toString() );
446
    if ( layerRef.at( 0 ) == QLatin1String( "\"" ) && layerRef.at( layerRef.size() - 1 ) == QLatin1String( "\"" ) )
447
    {
448
      layerRef.remove( 0, 1 );
449
      layerRef.chop( 1 );
450

    
451
    }
452
    layerRef.remove( QChar( '\\' ), Qt::CaseInsensitive );
453
    rasterReferences << layerRef;
454
  }
455

    
456
  return rasterReferences;
457
}
458

    
459
QgsRasterMatrix QgsRasterCalcNode::evaluateFunction( const QVector<QgsRasterMatrix *> &matrixVector, QgsRasterMatrix &result ) const
460
{
461

    
462
  if ( mFunctionName == "if" )
463
  {
464
    //scalar condition
465
    if ( matrixVector.at( 0 )->isNumber() )
466
    {
467
      result = ( matrixVector.at( 0 )->data() ? * matrixVector.at( 1 ) : * matrixVector.at( 2 ) );
468
      return result;
469
    }
470
    int nCols = matrixVector.at( 0 )->nColumns();
471
    int nRows = matrixVector.at( 0 )->nRows();
472
    int nEntries = nCols * nRows;
473
    std::unique_ptr< double > dataResult( new double[nEntries] );
474
    double *dataResultRawPtr =  dataResult.get();
475

    
476
    double *condition = matrixVector.at( 0 )->data();
477
    double *firstOption = matrixVector.at( 1 )->data();
478
    double *secondOption = matrixVector.at( 2 )->data();
479

    
480
    bool isFirstOptionNumber = matrixVector.at( 1 )->isNumber();
481
    bool isSecondCOptionNumber = matrixVector.at( 2 )->isNumber();
482
    double noDataValueCondition = matrixVector.at( 0 )->nodataValue();
483

    
484
    for ( int i = 0; i < nEntries; ++i )
485
    {
486
      if ( condition[i] == noDataValueCondition )
487
      {
488
        dataResultRawPtr[i] = result.nodataValue();
489
        continue;
490
      }
491
      else if ( condition[i] != 0 )
492
      {
493
        dataResultRawPtr[i] = isFirstOptionNumber ? firstOption[0] : firstOption[i];
494
        continue;
495
      }
496
      dataResultRawPtr[i] = isSecondCOptionNumber ? secondOption[0] : secondOption[i];
497
    }
498

    
499
    result.setData( nCols, nRows, dataResult.release(), result.nodataValue() );
500
  }
501
  return result;
502
}