Add if() function to raster calculator (#44839)

* start to work on new branch for conditional statement

* adjust the layout

* some pseudocode

* start to design the type tFunct, that should lead to the conditional statement

* modify the raw with a smart pointer

* change the test method and some other parts in the tFunct type

* complete the conditional statement option and update the test

* change evaluation method

* some optimization in the evaluation method

* minor adjustment

* minor adjustmentto test method

* add the button to the ui and some change to the code

* add a comment

* modify the parser and lexer in order to let the raster calc work with case-insensitive IF/if/If/iF

* change some parts according to the review and simplify the test method

* minor changes

* modify comment

* minor changes to enum type (tFunction)

* add some parts to test toString()  method

* add the possibility to use scalar condition in eveluationFunction() method and the corresponding test code

* update toString method

* update and optimize toString method

Co-authored-by: franc <Franc-Brs>
This commit is contained in:
Francesco Bursi 2021-09-07 09:48:12 +02:00 committed by GitHub
parent 5dc41385e2
commit d602f77a33
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
9 changed files with 291 additions and 3 deletions

View File

@ -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 <QgsRasterCalcNode *> functionArgs );
%Docstring
Constructor for the tFunction type
%End
QgsRasterCalcNode( const QString &rasterName );
~QgsRasterCalcNode();

View File

@ -73,6 +73,8 @@ raster_band_ref_quoted \"(\\.|[^"])*\"
"<=" { return LE; }
">=" { return GE; }
"if" { return IF; }
[=><+-/*^] { return yytext[0]; }

View File

@ -36,6 +36,13 @@ QgsRasterCalcNode::QgsRasterCalcNode( Operator op, QgsRasterCalcNode *left, QgsR
{
}
QgsRasterCalcNode::QgsRasterCalcNode( QString functionName, QVector <QgsRasterCalcNode *> functionArgs )
: mType( tFunction )
, mFunctionName( functionName )
, mFunctionArgs( functionArgs )
{
}
QgsRasterCalcNode::QgsRasterCalcNode( const QString &rasterName )
: mType( tRasterRef )
, mRasterName( rasterName )
@ -208,6 +215,21 @@ bool QgsRasterCalcNode::calculate( QMap<QString, QgsRasterBlock * > &rasterData,
result.setData( mMatrix->nColumns(), mMatrix->nRows(), data, result.nodataValue() );
return true;
}
else if ( mType == tFunction )
{
QVector <QgsRasterMatrix *> 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<const QgsRasterCalcNode *> 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<QgsRasterMatrix *> &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;
}

View File

@ -25,6 +25,8 @@
#include <QString>
#include "qgis_analysis.h"
#include <QVector>
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 <QgsRasterCalcNode *> 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<QgsRasterMatrix *> &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 <QgsRasterCalcNode *> mFunctionArgs;
};

View File

@ -58,6 +58,8 @@
%token<op> FUNCTION
%token<op> FUNCTION_2_ARGS
%token IF
%type <node> root
%type <node> 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 <QgsRasterCalcNode *> 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); }

View File

@ -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 )
{
// '"' -> '\\"'

View File

@ -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

View File

@ -492,6 +492,13 @@
</property>
</widget>
</item>
<item row="3" column="4">
<widget class="QPushButton" name="mConditionalStatButton">
<property name="text">
<string notr="true">if</string>
</property>
</widget>
</item>
</layout>
</widget>
</item>

View File

@ -26,6 +26,8 @@ Email : nyall dot dawson at gmail dot com
#include "qgsapplication.h"
#include "qgsproject.h"
#include <QDebug>
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<QgsRasterCalculatorEntry> 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"