diff --git a/python/analysis/auto_generated/raster/qgsrastercalcnode.sip.in b/python/analysis/auto_generated/raster/qgsrastercalcnode.sip.in index b073178a941..9f3dacca2b8 100644 --- a/python/analysis/auto_generated/raster/qgsrastercalcnode.sip.in +++ b/python/analysis/auto_generated/raster/qgsrastercalcnode.sip.in @@ -10,6 +10,7 @@ + class QgsRasterCalcNode { %Docstring(signature="appended") @@ -25,7 +26,8 @@ Represents a node in a raster calculator. tOperator, tNumber, tRasterRef, - tMatrix + tMatrix, + tFunction }; enum Operator @@ -67,6 +69,10 @@ Constructor for QgsRasterCalcNode. QgsRasterCalcNode( double number ); QgsRasterCalcNode( QgsRasterMatrix *matrix ); QgsRasterCalcNode( Operator op, QgsRasterCalcNode *left, QgsRasterCalcNode *right ); + QgsRasterCalcNode( QString functionName, QVector functionArgs ); +%Docstring +Constructor for the tFunction type +%End QgsRasterCalcNode( const QString &rasterName ); ~QgsRasterCalcNode(); diff --git a/src/analysis/raster/qgsrastercalclexer.ll b/src/analysis/raster/qgsrastercalclexer.ll index a10caf154e5..fb074f4016c 100644 --- a/src/analysis/raster/qgsrastercalclexer.ll +++ b/src/analysis/raster/qgsrastercalclexer.ll @@ -73,6 +73,8 @@ raster_band_ref_quoted \"(\\.|[^"])*\" "<=" { return LE; } ">=" { return GE; } +"if" { return IF; } + [=><+-/*^] { return yytext[0]; } diff --git a/src/analysis/raster/qgsrastercalcnode.cpp b/src/analysis/raster/qgsrastercalcnode.cpp index 4a653cf8efa..f8328e78f39 100644 --- a/src/analysis/raster/qgsrastercalcnode.cpp +++ b/src/analysis/raster/qgsrastercalcnode.cpp @@ -36,6 +36,13 @@ QgsRasterCalcNode::QgsRasterCalcNode( Operator op, QgsRasterCalcNode *left, QgsR { } +QgsRasterCalcNode::QgsRasterCalcNode( QString functionName, QVector functionArgs ) + : mType( tFunction ) + , mFunctionName( functionName ) + , mFunctionArgs( functionArgs ) +{ +} + QgsRasterCalcNode::QgsRasterCalcNode( const QString &rasterName ) : mType( tRasterRef ) , mRasterName( rasterName ) @@ -208,6 +215,21 @@ bool QgsRasterCalcNode::calculate( QMap &rasterData, result.setData( mMatrix->nColumns(), mMatrix->nRows(), data, result.nodataValue() ); return true; } + else if ( mType == tFunction ) + { + QVector matrixContainer; + for ( int i = 0; i < mFunctionArgs.size(); ++i ) + { + std::unique_ptr< QgsRasterMatrix > singleMatrix( new QgsRasterMatrix( result.nColumns(), result.nRows(), nullptr, result.nodataValue() ) ); + if ( !mFunctionArgs.at( i ) || !mFunctionArgs.at( i )->calculate( rasterData, *singleMatrix, row ) ) + { + return false; + } + matrixContainer.append( singleMatrix.release() ); + } + evaluateFunction( matrixContainer, result ); + return true; + } return false; } @@ -360,6 +382,18 @@ QString QgsRasterCalcNode::toString( bool cStyle ) const break; case tMatrix: break; + case tFunction: + if ( mFunctionName == "if" ) + { + const QString argOne = mFunctionArgs.at( 0 )->toString( cStyle ); + const QString argTwo = mFunctionArgs.at( 1 )->toString( cStyle ); + const QString argThree = mFunctionArgs.at( 2 )->toString( cStyle ); + if ( cStyle ) + result = QStringLiteral( " ( %1 ) ? ( %2 ) : ( %3 ) " ).arg( argOne, argTwo, argThree ); + else + result = QStringLiteral( "if( %1 , %2 , %3 )" ).arg( argOne, argTwo, argThree ); + } + break; } return result; } @@ -373,6 +407,10 @@ QList QgsRasterCalcNode::findNodes( const QgsRasterCa nodeList.append( mLeft->findNodes( type ) ); if ( mRight ) nodeList.append( mRight->findNodes( type ) ); + + for ( QgsRasterCalcNode *node : mFunctionArgs ) + nodeList.append( node->findNodes( type ) ); + return nodeList; } @@ -417,3 +455,48 @@ QStringList QgsRasterCalcNode::cleanRasterReferences() return rasterReferences; } + +QgsRasterMatrix QgsRasterCalcNode::evaluateFunction( const QVector &matrixVector, QgsRasterMatrix &result ) const +{ + + if ( mFunctionName == "if" ) + { + //scalar condition + if ( matrixVector.at( 0 )->isNumber() ) + { + result = ( matrixVector.at( 0 )->data() ? * matrixVector.at( 1 ) : * matrixVector.at( 2 ) ); + return result; + } + int nCols = matrixVector.at( 0 )->nColumns(); + int nRows = matrixVector.at( 0 )->nRows(); + int nEntries = nCols * nRows; + std::unique_ptr< double > dataResult( new double[nEntries] ); + double *dataResultRawPtr = dataResult.get(); + + double *condition = matrixVector.at( 0 )->data(); + double *firstOption = matrixVector.at( 1 )->data(); + double *secondOption = matrixVector.at( 2 )->data(); + + bool isFirstOptionNumber = matrixVector.at( 1 )->isNumber(); + bool isSecondCOptionNumber = matrixVector.at( 2 )->isNumber(); + double noDataValueCondition = matrixVector.at( 0 )->nodataValue(); + + for ( int i = 0; i < nEntries; ++i ) + { + if ( condition[i] == noDataValueCondition ) + { + dataResultRawPtr[i] = result.nodataValue(); + continue; + } + else if ( condition[i] != 0 ) + { + dataResultRawPtr[i] = isFirstOptionNumber ? firstOption[0] : firstOption[i]; + continue; + } + dataResultRawPtr[i] = isSecondCOptionNumber ? secondOption[0] : secondOption[i]; + } + + result.setData( nCols, nRows, dataResult.release(), result.nodataValue() ); + } + return result; +} diff --git a/src/analysis/raster/qgsrastercalcnode.h b/src/analysis/raster/qgsrastercalcnode.h index 8b20cca088d..331eeefcf63 100644 --- a/src/analysis/raster/qgsrastercalcnode.h +++ b/src/analysis/raster/qgsrastercalcnode.h @@ -25,6 +25,8 @@ #include #include "qgis_analysis.h" +#include + class QgsRasterBlock; class QgsRasterMatrix; @@ -42,7 +44,8 @@ class ANALYSIS_EXPORT QgsRasterCalcNode tOperator = 1, tNumber, tRasterRef, - tMatrix + tMatrix, + tFunction }; //! possible operators @@ -85,6 +88,8 @@ class ANALYSIS_EXPORT QgsRasterCalcNode QgsRasterCalcNode( double number ); QgsRasterCalcNode( QgsRasterMatrix *matrix ); QgsRasterCalcNode( Operator op, QgsRasterCalcNode *left, QgsRasterCalcNode *right ); + //!Constructor for the tFunction type + QgsRasterCalcNode( QString functionName, QVector functionArgs ); QgsRasterCalcNode( const QString &rasterName ); ~QgsRasterCalcNode(); @@ -143,6 +148,12 @@ class ANALYSIS_EXPORT QgsRasterCalcNode QgsRasterCalcNode( const QgsRasterCalcNode &rh ); #endif + /** + * Calculates result of raster calculation when tFunct type is used + * \since QGIS 3.22 + */ + QgsRasterMatrix evaluateFunction( const QVector &matrixVector, QgsRasterMatrix &result ) const; + Type mType = tNumber; QgsRasterCalcNode *mLeft = nullptr; QgsRasterCalcNode *mRight = nullptr; @@ -150,7 +161,9 @@ class ANALYSIS_EXPORT QgsRasterCalcNode QString mRasterName; QgsRasterMatrix *mMatrix = nullptr; Operator mOperator = opNONE; - + //added for the conditional statement + QString mFunctionName; + QVector mFunctionArgs; }; diff --git a/src/analysis/raster/qgsrastercalcparser.yy b/src/analysis/raster/qgsrastercalcparser.yy index 8b40b61191b..dfd1d3321ca 100644 --- a/src/analysis/raster/qgsrastercalcparser.yy +++ b/src/analysis/raster/qgsrastercalcparser.yy @@ -58,6 +58,8 @@ %token FUNCTION %token FUNCTION_2_ARGS +%token IF + %type root %type raster_exp @@ -81,6 +83,12 @@ root: raster_exp{} raster_exp: FUNCTION '(' raster_exp ')' { $$ = new QgsRasterCalcNode($1, $3, 0); joinTmpNodes($$, $3, 0);} | FUNCTION_2_ARGS '(' raster_exp ',' raster_exp ')' { $$ = new QgsRasterCalcNode($1, $3, $5); joinTmpNodes($$, $3, $5);} + | IF '(' raster_exp ',' raster_exp ',' raster_exp ')' { QVector tmpVect; + tmpVect<< $3<< $5<< $7; + $$ = new QgsRasterCalcNode("if", tmpVect); + joinTmpNodes($$, $3, $5); + gTmpNodes.removeAll($7); + } | raster_exp AND raster_exp { $$ = new QgsRasterCalcNode( QgsRasterCalcNode::opAND, $1, $3 ); joinTmpNodes($$,$1,$3); } | raster_exp OR raster_exp { $$ = new QgsRasterCalcNode( QgsRasterCalcNode::opOR, $1, $3 ); joinTmpNodes($$,$1,$3); } | raster_exp '=' raster_exp { $$ = new QgsRasterCalcNode( QgsRasterCalcNode::opEQ, $1, $3 ); joinTmpNodes($$,$1,$3); } diff --git a/src/app/qgsrastercalcdialog.cpp b/src/app/qgsrastercalcdialog.cpp index b8cbc0a3257..07ad718a71c 100644 --- a/src/app/qgsrastercalcdialog.cpp +++ b/src/app/qgsrastercalcdialog.cpp @@ -73,6 +73,7 @@ QgsRasterCalcDialog::QgsRasterCalcDialog( QgsRasterLayer *rasterLayer, QWidget * connect( mMinButton, &QPushButton::clicked, this, &QgsRasterCalcDialog::mMinButton_clicked ); connect( mMaxButton, &QPushButton::clicked, this, &QgsRasterCalcDialog::mMaxButton_clicked ); connect( mOrButton, &QPushButton::clicked, this, &QgsRasterCalcDialog::mOrButton_clicked ); + connect( mConditionalStatButton, &QPushButton::clicked, this, &QgsRasterCalcDialog::mConditionalStatButton_clicked ); connect( mButtonBox, &QDialogButtonBox::helpRequested, this, &QgsRasterCalcDialog::showHelp ); if ( rasterLayer && rasterLayer->dataProvider() && rasterLayer->providerType() == QLatin1String( "gdal" ) ) @@ -539,6 +540,11 @@ void QgsRasterCalcDialog::mMaxButton_clicked() mExpressionTextEdit->insertPlainText( QStringLiteral( " MAX ( " ) ); } +void QgsRasterCalcDialog::mConditionalStatButton_clicked() +{ + mExpressionTextEdit->insertPlainText( QStringLiteral( " if ( " ) ); +} + QString QgsRasterCalcDialog::quoteBandEntry( const QString &layerName ) { // '"' -> '\\"' diff --git a/src/app/qgsrastercalcdialog.h b/src/app/qgsrastercalcdialog.h index 7904c0b6cbc..590806ba7d1 100644 --- a/src/app/qgsrastercalcdialog.h +++ b/src/app/qgsrastercalcdialog.h @@ -100,6 +100,7 @@ class APP_EXPORT QgsRasterCalcDialog: public QDialog, private Ui::QgsRasterCalcD void mAbsButton_clicked(); void mMinButton_clicked(); void mMaxButton_clicked(); + void mConditionalStatButton_clicked(); private: //! Sets the extent and size of the output diff --git a/src/ui/qgsrastercalcdialogbase.ui b/src/ui/qgsrastercalcdialogbase.ui index ee37a51665b..3245b32d51d 100644 --- a/src/ui/qgsrastercalcdialogbase.ui +++ b/src/ui/qgsrastercalcdialogbase.ui @@ -492,6 +492,13 @@ + + + + if + + + diff --git a/tests/src/analysis/testqgsrastercalculator.cpp b/tests/src/analysis/testqgsrastercalculator.cpp index dfb35c59412..3afd8bf2c68 100644 --- a/tests/src/analysis/testqgsrastercalculator.cpp +++ b/tests/src/analysis/testqgsrastercalculator.cpp @@ -26,6 +26,8 @@ Email : nyall dot dawson at gmail dot com #include "qgsapplication.h" #include "qgsproject.h" +#include + Q_DECLARE_METATYPE( QgsRasterCalcNode::Operator ) class TestQgsRasterCalculator : public QObject @@ -67,6 +69,9 @@ class TestQgsRasterCalculator : public QObject void testStatistics(); + void parseFunctionTypeString(); //test the parsing of the formule for the tFunction type + void testFunctionTypeWithLayer(); //test of conditional statement + private: QgsRasterLayer *mpLandsatRasterLayer = nullptr; @@ -744,6 +749,9 @@ void TestQgsRasterCalculator::toString() // Test regression #32477 QCOMPARE( _test( QStringLiteral( R"raw(("r@1"<100.09)*0.1)raw" ), true ), QString( R"raw(( float ) ( ( float ) "r@1" < ( float ) 100.09 ) * ( float ) 0.1)raw" ) ); + //test the conditional statement + QCOMPARE( _test( QStringLiteral( "if( \"raster@1\" > 5 , 100 , 5)" ), false ), QString( "if( \"raster@1\" > 5 , 100 , 5 )" ) ); + QCOMPARE( _test( QStringLiteral( "if( \"raster@1\" > 5 , 100 , 5)" ), true ), QString( " ( ( float ) ( ( float ) \"raster@1\" > ( float ) 5 ) ) ? ( ( float ) 100 ) : ( ( float ) 5 ) " ) ); QString error; std::unique_ptr< QgsRasterCalcNode > calcNode( QgsRasterCalcNode::parseRasterCalcString( QStringLiteral( "min( \"raster@1\" )" ), error ) ); @@ -876,6 +884,160 @@ void TestQgsRasterCalculator::testStatistics() } +void TestQgsRasterCalculator::parseFunctionTypeString() +{ + QString errorString; + const QgsRasterCalcNode *node { QgsRasterCalcNode::parseRasterCalcString( QString( ), errorString ) }; + QVERIFY( ! node ); + QVERIFY( ! errorString.isEmpty() ); + + errorString = QString(); + node = QgsRasterCalcNode::parseRasterCalcString( QStringLiteral( "if(\"raster@1\">5,100,5)" ), errorString ); + QVERIFY( node ); + QVERIFY( errorString.isEmpty() ); + QVERIFY( node->findNodes( QgsRasterCalcNode::Type::tRasterRef ).length() == 1 ); + + //test case sensitivity (instead of "if", use "IF") + errorString = QString(); + node = QgsRasterCalcNode::parseRasterCalcString( QStringLiteral( "IF(\"raster@1\">5,100,5)" ), errorString ); + QVERIFY( node ); + QVERIFY( errorString.isEmpty() ); + QVERIFY( node->findNodes( QgsRasterCalcNode::Type::tRasterRef ).length() == 1 ); +} + +void TestQgsRasterCalculator::testFunctionTypeWithLayer() +{ + // first band + QgsRasterCalculatorEntry entry1; + entry1.bandNumber = 1; + entry1.raster = mpLandsatRasterLayer; + entry1.ref = QStringLiteral( "landsat@1" ); + + // second band + QgsRasterCalculatorEntry entry2; + entry2.bandNumber = 2; + entry2.raster = mpLandsatRasterLayer; + entry2.ref = QStringLiteral( "landsat@2" ); + + QVector entries; + entries << entry1 << entry2; + + QgsCoordinateReferenceSystem crs( QStringLiteral( "EPSG:32633" ) ); + QgsRectangle extent( 783235, 3348110, 783350, 3347960 ); + + QTemporaryFile tmpFile; + tmpFile.open(); // fileName is not available until open + QString tmpName = tmpFile.fileName(); + tmpFile.close(); + + // Test with one raster as condition and numbers as first and second option + // if ( landsat@1 > 124.5, 100.0 , 5.0 ) + QgsRasterCalculator rc( QStringLiteral( " if(\"landsat@1\">124.5, 100.0 , 5.0 ) " ), + tmpName, + QStringLiteral( "GTiff" ), + extent, crs, 2, 3, entries, + QgsProject::instance()->transformContext() ); + QCOMPARE( static_cast< int >( rc.processCalculation() ), 0 ); + + //open output file and check results + QgsRasterLayer *result = new QgsRasterLayer( tmpName, QStringLiteral( "result" ) ); + QCOMPARE( result->width(), 2 ); + QCOMPARE( result->height(), 3 ); + QgsRasterBlock *block = result->dataProvider()->block( 1, extent, 2, 3 ); + + QCOMPARE( block->value( 0, 0 ), 100.0 ); + QCOMPARE( block->value( 0, 1 ), 100.0 ); + QCOMPARE( block->value( 1, 0 ), 5.0 ); + QCOMPARE( block->value( 1, 1 ), 100.0 ); + QCOMPARE( block->value( 2, 0 ), 100.0 ); + QCOMPARE( block->value( 2, 1 ), 5.0 ); + + // Test with one raster as condition, one raster first option and number as second option + // if ( landsat@1 > 124.5, landsat@1 + landsat@2 , 5.0 ) + QgsRasterCalculator rc2( QStringLiteral( " if(\"landsat@1\">124.5, \"landsat@1\" + \"landsat@2\" , 5.0 ) " ), + tmpName, + QStringLiteral( "GTiff" ), + extent, crs, 2, 3, entries, + QgsProject::instance()->transformContext() ); + QCOMPARE( static_cast< int >( rc2.processCalculation() ), 0 ); + + //open output file and check results + result = new QgsRasterLayer( tmpName, QStringLiteral( "result" ) ); + QCOMPARE( result->width(), 2 ); + QCOMPARE( result->height(), 3 ); + block = result->dataProvider()->block( 1, extent, 2, 3 ); + QCOMPARE( block->value( 0, 0 ), 265.0 ); + QCOMPARE( block->value( 0, 1 ), 263.0 ); + QCOMPARE( block->value( 1, 0 ), 5.0 ); + QCOMPARE( block->value( 1, 1 ), 264.0 ); + QCOMPARE( block->value( 2, 0 ), 266.0 ); + QCOMPARE( block->value( 2, 1 ), 5.0 ); + + // Test with one raster as condition, one raster first option and number as second option + // if ( landsat@1 > 124.5, landsat@1 + landsat@2 , landsat@3 ) + QgsRasterCalculator rc3( QStringLiteral( " if(\"landsat@1\">124.5, \"landsat@1\" + \"landsat@2\" , \"landsat@1\" - \"landsat@2\" ) " ), + tmpName, + QStringLiteral( "GTiff" ), + extent, crs, 2, 3, entries, + QgsProject::instance()->transformContext() ); + QCOMPARE( static_cast< int >( rc3.processCalculation() ), 0 ); + + //open output file and check results + result = new QgsRasterLayer( tmpName, QStringLiteral( "result" ) ); + QCOMPARE( result->width(), 2 ); + QCOMPARE( result->height(), 3 ); + block = result->dataProvider()->block( 1, extent, 2, 3 ); + QCOMPARE( block->value( 0, 0 ), 265.0 ); + QCOMPARE( block->value( 0, 1 ), 263.0 ); + QCOMPARE( block->value( 1, 0 ), -15.0 ); + QCOMPARE( block->value( 1, 1 ), 264.0 ); + QCOMPARE( block->value( 2, 0 ), 266.0 ); + QCOMPARE( block->value( 2, 1 ), -13.0 ); + + // Test with scalar (always true) as condition, one raster first option and number as second option + // if ( 5 > 4, landsat@1 + landsat@2 , 0 ) + QgsRasterCalculator rc4( QStringLiteral( " if( 5>4 , \"landsat@1\" + \"landsat@2\" , 0 ) " ), + tmpName, + QStringLiteral( "GTiff" ), + extent, crs, 2, 3, entries, + QgsProject::instance()->transformContext() ); + QCOMPARE( static_cast< int >( rc4.processCalculation() ), 0 ); + + //open output file and check results + result = new QgsRasterLayer( tmpName, QStringLiteral( "result" ) ); + QCOMPARE( result->width(), 2 ); + QCOMPARE( result->height(), 3 ); + block = result->dataProvider()->block( 1, extent, 2, 3 ); + QCOMPARE( block->value( 0, 0 ), 265.0 ); + QCOMPARE( block->value( 0, 1 ), 263.0 ); + QCOMPARE( block->value( 1, 0 ), 263.0 ); + QCOMPARE( block->value( 1, 1 ), 264.0 ); + QCOMPARE( block->value( 2, 0 ), 266.0 ); + QCOMPARE( block->value( 2, 1 ), 261.0 ); + + // Test with scalar (always false) as condition, one raster first option and number as second option + // if ( 4 > 5, landsat@1 + landsat@2 , 0 ) + QgsRasterCalculator rc5( QStringLiteral( " if( 4>5 , \"landsat@1\" + \"landsat@2\" , 0 ) " ), + tmpName, + QStringLiteral( "GTiff" ), + extent, crs, 2, 3, entries, + QgsProject::instance()->transformContext() ); + QCOMPARE( static_cast< int >( rc5.processCalculation() ), 0 ); + + //open output file and check results + result = new QgsRasterLayer( tmpName, QStringLiteral( "result" ) ); + QCOMPARE( result->width(), 2 ); + QCOMPARE( result->height(), 3 ); + block = result->dataProvider()->block( 1, extent, 2, 3 ); + QCOMPARE( block->value( 0, 0 ), 0.0 ); + QCOMPARE( block->value( 0, 1 ), 0.0 ); + QCOMPARE( block->value( 1, 0 ), 0.0 ); + QCOMPARE( block->value( 1, 1 ), 0.0 ); + QCOMPARE( block->value( 2, 0 ), 0.0 ); + QCOMPARE( block->value( 2, 1 ), 0.0 ); + delete result; + delete block; +} QGSTEST_MAIN( TestQgsRasterCalculator ) #include "testqgsrastercalculator.moc"