Skip to content

Commit a0e2ee4

Browse files
elpasoelpaso
authoredJun 15, 2021
Merge pull request #43512 from elpaso/bugfix-gh42835-rastercalc-missing-stats
Fix missing stats after raster calc output
2 parents 51c0886 + c1fa7c9 commit a0e2ee4

File tree

3 files changed

+87
-14
lines changed

3 files changed

+87
-14
lines changed
 

‎src/analysis/raster/qgsrastercalculator.cpp

Lines changed: 45 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,41 @@
3636
#include "qgsgdalutils.h"
3737
#endif
3838

39+
40+
//
41+
// global callback function
42+
//
43+
int CPL_STDCALL GdalProgressCallback( double dfComplete,
44+
const char *pszMessage,
45+
void *pProgressArg )
46+
{
47+
Q_UNUSED( pszMessage )
48+
49+
static double sDfLastComplete = -1.0;
50+
51+
QgsFeedback *feedback = static_cast<QgsFeedback *>( pProgressArg );
52+
53+
if ( sDfLastComplete > dfComplete )
54+
{
55+
if ( sDfLastComplete >= 1.0 )
56+
sDfLastComplete = -1.0;
57+
else
58+
sDfLastComplete = dfComplete;
59+
}
60+
61+
if ( std::floor( sDfLastComplete * 10 ) != std::floor( dfComplete * 10 ) )
62+
{
63+
if ( feedback )
64+
feedback->setProgress( dfComplete * 100 );
65+
}
66+
sDfLastComplete = dfComplete;
67+
68+
if ( feedback && feedback->isCanceled() )
69+
return false;
70+
71+
return true;
72+
}
73+
3974
QgsRasterCalculator::QgsRasterCalculator( const QString &formulaString, const QString &outputFile, const QString &outputFormat, const QgsRectangle &outputExtent, int nOutputColumns, int nOutputRows, const QVector<QgsRasterCalculatorEntry> &rasterEntries, const QgsCoordinateTransformContext &transformContext )
4075
: mFormulaString( formulaString )
4176
, mOutputFile( outputFile )
@@ -131,7 +166,6 @@ QgsRasterCalculator::Result QgsRasterCalculator::processCalculation( QgsFeedback
131166

132167
#ifdef HAVE_OPENCL
133168
// Check for matrix nodes, GPU implementation does not support them
134-
QList<const QgsRasterCalcNode *> nodeList;
135169
if ( QgsOpenClUtils::enabled() && QgsOpenClUtils::available() && ! requiresMatrix )
136170
{
137171
return processCalculationGPU( std::move( calcNode ), feedback );
@@ -354,6 +388,9 @@ QgsRasterCalculator::Result QgsRasterCalculator::processCalculation( QgsFeedback
354388
gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
355389
return Canceled;
356390
}
391+
392+
GDALComputeRasterStatistics( outputRasterBand, true, nullptr, nullptr, nullptr, nullptr, GdalProgressCallback, feedback );
393+
357394
return Success;
358395
}
359396

@@ -399,7 +436,7 @@ QgsRasterCalculator::Result QgsRasterCalculator::processCalculationGPU( std::uni
399436
LayerRef entry;
400437
entry.name = r;
401438
entry.band = parts[1].toInt( &ok );
402-
for ( const auto &ref : mRasterEntries )
439+
for ( const auto &ref : std::as_const( mRasterEntries ) )
403440
{
404441
if ( ref.ref == entry.name )
405442
entry.layer = ref.raster;
@@ -459,8 +496,7 @@ QgsRasterCalculator::Result QgsRasterCalculator::processCalculationGPU( std::uni
459496
{
460497
cExpression.replace( QStringLiteral( "\"%1\"" ).arg( ref.name ), QStringLiteral( "%1[i]" ).arg( ref.varName ) );
461498
inputArgs.append( QStringLiteral( "__global %1 *%2" )
462-
.arg( ref.typeName )
463-
.arg( ref.varName ) );
499+
.arg( ref.typeName, ref.varName ) );
464500
inputBuffers.push_back( cl::Buffer( ctx, CL_MEM_READ_ONLY, ref.bufferSize, nullptr, nullptr ) );
465501
}
466502

@@ -490,10 +526,10 @@ QgsRasterCalculator::Result QgsRasterCalculator::processCalculationGPU( std::uni
490526
QStringList inputNoDataCheck;
491527
for ( const auto &ref : inputRefs )
492528
{
493-
inputDesc.append( QStringLiteral( " // %1 = %2" ).arg( ref.varName ).arg( ref.name ) );
529+
inputDesc.append( QStringLiteral( " // %1 = %2" ).arg( ref.varName, ref.name ) );
494530
if ( ref.layer->dataProvider()->sourceHasNoDataValue( ref.band ) )
495531
{
496-
inputNoDataCheck.append( QStringLiteral( "(float) %1[i] == (float) %2" ).arg( ref.varName ).arg( ref.layer->dataProvider()->sourceNoDataValue( ref.band ) ) );
532+
inputNoDataCheck.append( QStringLiteral( "(float) %1[i] == (float) %2" ).arg( ref.varName, QString::number( ref.layer->dataProvider()->sourceNoDataValue( ref.band ), 'g', 10 ) ) );
497533
}
498534
}
499535

@@ -503,7 +539,7 @@ QgsRasterCalculator::Result QgsRasterCalculator::processCalculationGPU( std::uni
503539
programTemplate = programTemplate.replace( QLatin1String( "##EXPRESSION##" ), cExpression );
504540
programTemplate = programTemplate.replace( QLatin1String( "##EXPRESSION_ORIGINAL##" ), calcNode->toString( ) );
505541

506-
// qDebug() << programTemplate;
542+
//qDebug() << programTemplate;
507543

508544
// Create a program from the kernel source
509545
cl::Program program( QgsOpenClUtils::buildProgram( programTemplate, QgsOpenClUtils::ExceptionBehavior::Throw ) );
@@ -628,6 +664,8 @@ QgsRasterCalculator::Result QgsRasterCalculator::processCalculationGPU( std::uni
628664

629665
inputBuffers.clear();
630666

667+
GDALComputeRasterStatistics( outputRasterBand, true, nullptr, nullptr, nullptr, nullptr, GdalProgressCallback, feedback );
668+
631669
}
632670
catch ( cl::Error &e )
633671
{

‎src/app/qgisapp.cpp

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6816,6 +6816,7 @@ void QgisApp::showRasterCalculator()
68166816
QgsProject::instance()->transformContext() );
68176817

68186818
QProgressDialog p( tr( "Calculating raster expression…" ), tr( "Abort" ), 0, 0 );
6819+
p.setWindowTitle( tr( "Raster calculator" ) );
68196820
p.setWindowModality( Qt::WindowModal );
68206821
p.setMaximum( 100.0 );
68216822
QgsFeedback feedback;

‎tests/src/analysis/testqgsrastercalculator.cpp

Lines changed: 41 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -65,6 +65,8 @@ class TestQgsRasterCalculator : public QObject
6565
void testRasterEntries();
6666
void calcFormulasWithReprojectedLayers();
6767

68+
void testStatistics();
69+
6870
private:
6971

7072
QgsRasterLayer *mpLandsatRasterLayer = nullptr;
@@ -799,13 +801,13 @@ void TestQgsRasterCalculator::calcFormulasWithReprojectedLayers()
799801
QgsRasterBlock *block = result->dataProvider()->block( 1, extent, 2, 3 );
800802
qDebug() << "Actual:" << block->value( 0, 0 ) << block->value( 0, 1 ) << block->value( 1, 0 ) << block->value( 1, 1 ) << block->value( 2, 0 ) << block->value( 2, 1 );
801803
qDebug() << "Expected:" << values[0] << values[1] << values[2] << values[3] << values[4] << values[5];
802-
const double epsilon { 0.0001 };
803-
QVERIFY( qgsDoubleNear( block->value( 0, 0 ), static_cast<double>( values[0] ), epsilon ) );
804-
QVERIFY( qgsDoubleNear( block->value( 0, 1 ), static_cast<double>( values[1] ), epsilon ) );
805-
QVERIFY( qgsDoubleNear( block->value( 1, 0 ), static_cast<double>( values[2] ), epsilon ) );
806-
QVERIFY( qgsDoubleNear( block->value( 1, 1 ), static_cast<double>( values[3] ), epsilon ) );
807-
QVERIFY( qgsDoubleNear( block->value( 2, 0 ), static_cast<double>( values[4] ), epsilon ) );
808-
QVERIFY( qgsDoubleNear( block->value( 2, 1 ), static_cast<double>( values[5] ), epsilon ) );
804+
const double epsilon { 0.001 };
805+
QVERIFY2( qgsDoubleNear( block->value( 0, 0 ), static_cast<double>( values[0] ), epsilon ), formula.toUtf8().constData() );
806+
QVERIFY2( qgsDoubleNear( block->value( 0, 1 ), static_cast<double>( values[1] ), epsilon ), formula.toUtf8().constData() );
807+
QVERIFY2( qgsDoubleNear( block->value( 1, 0 ), static_cast<double>( values[2] ), epsilon ), formula.toUtf8().constData() );
808+
QVERIFY2( qgsDoubleNear( block->value( 1, 1 ), static_cast<double>( values[3] ), epsilon ), formula.toUtf8().constData() );
809+
QVERIFY2( qgsDoubleNear( block->value( 2, 0 ), static_cast<double>( values[4] ), epsilon ), formula.toUtf8().constData() );
810+
QVERIFY2( qgsDoubleNear( block->value( 2, 1 ), static_cast<double>( values[5] ), epsilon ), formula.toUtf8().constData() );
809811
delete result;
810812
delete block;
811813
};
@@ -842,6 +844,38 @@ void TestQgsRasterCalculator::calcFormulasWithReprojectedLayers()
842844

843845
}
844846

847+
void TestQgsRasterCalculator::testStatistics()
848+
{
849+
QgsRasterCalculatorEntry entry1;
850+
entry1.bandNumber = 1;
851+
entry1.raster = mpLandsatRasterLayer;
852+
entry1.ref = QStringLiteral( "landsat@1" );
853+
854+
QTemporaryFile tmpFile;
855+
tmpFile.open(); // fileName is not available until open
856+
QString tmpName = tmpFile.fileName();
857+
tmpFile.close();
858+
859+
QgsCoordinateReferenceSystem crs( QStringLiteral( "EPSG:32633" ) );
860+
QgsRectangle extent( 783235, 3348110, 783350, 3347960 );
861+
862+
QgsRasterCalculator rc( QStringLiteral( "\"landsat@1\" * 2" ),
863+
tmpName,
864+
QStringLiteral( "GTiff" ),
865+
extent, crs, 2, 3, { entry1 },
866+
QgsProject::instance()->transformContext() );
867+
QCOMPARE( static_cast< int >( rc.processCalculation() ), 0 );
868+
869+
//open output file and check stats are there
870+
auto ds = GDALOpenEx( tmpName.toUtf8().constData(), GDAL_OF_RASTER | GDAL_OF_READONLY, nullptr, nullptr, nullptr );
871+
auto band = GDALGetRasterBand( ds, 1 );
872+
double sMin, sMax, sMean, sStdDev;
873+
QCOMPARE( GDALGetRasterStatistics( band, true, false, &sMin, &sMax, &sMean, &sStdDev ), CE_None );
874+
QCOMPARE( sMin, 248.0 );
875+
QCOMPARE( sMax, 250.0 );
876+
877+
}
878+
845879

846880
QGSTEST_MAIN( TestQgsRasterCalculator )
847881
#include "testqgsrastercalculator.moc"

0 commit comments

Comments
 (0)
Please sign in to comment.