Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Fix missing stats after raster calc output
Fixes #42835

Also on raster calc:

- fix unreported issue with nodata float opencl
- papercut: set progress dialog window title
  • Loading branch information
elpaso committed Jun 2, 2021
1 parent c6fdec8 commit c1fa7c9
Show file tree
Hide file tree
Showing 3 changed files with 87 additions and 14 deletions.
52 changes: 45 additions & 7 deletions src/analysis/raster/qgsrastercalculator.cpp
Expand Up @@ -36,6 +36,41 @@
#include "qgsgdalutils.h"
#endif


//
// global callback function
//
int CPL_STDCALL GdalProgressCallback( double dfComplete,
const char *pszMessage,
void *pProgressArg )
{
Q_UNUSED( pszMessage )

static double sDfLastComplete = -1.0;

QgsFeedback *feedback = static_cast<QgsFeedback *>( pProgressArg );

if ( sDfLastComplete > dfComplete )
{
if ( sDfLastComplete >= 1.0 )
sDfLastComplete = -1.0;
else
sDfLastComplete = dfComplete;
}

if ( std::floor( sDfLastComplete * 10 ) != std::floor( dfComplete * 10 ) )
{
if ( feedback )
feedback->setProgress( dfComplete * 100 );
}
sDfLastComplete = dfComplete;

if ( feedback && feedback->isCanceled() )
return false;

return true;
}

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 )
: mFormulaString( formulaString )
, mOutputFile( outputFile )
Expand Down Expand Up @@ -131,7 +166,6 @@ QgsRasterCalculator::Result QgsRasterCalculator::processCalculation( QgsFeedback

#ifdef HAVE_OPENCL
// Check for matrix nodes, GPU implementation does not support them
QList<const QgsRasterCalcNode *> nodeList;
if ( QgsOpenClUtils::enabled() && QgsOpenClUtils::available() && ! requiresMatrix )
{
return processCalculationGPU( std::move( calcNode ), feedback );
Expand Down Expand Up @@ -354,6 +388,9 @@ QgsRasterCalculator::Result QgsRasterCalculator::processCalculation( QgsFeedback
gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
return Canceled;
}

GDALComputeRasterStatistics( outputRasterBand, true, nullptr, nullptr, nullptr, nullptr, GdalProgressCallback, feedback );

return Success;
}

Expand Down Expand Up @@ -399,7 +436,7 @@ QgsRasterCalculator::Result QgsRasterCalculator::processCalculationGPU( std::uni
LayerRef entry;
entry.name = r;
entry.band = parts[1].toInt( &ok );
for ( const auto &ref : mRasterEntries )
for ( const auto &ref : std::as_const( mRasterEntries ) )
{
if ( ref.ref == entry.name )
entry.layer = ref.raster;
Expand Down Expand Up @@ -459,8 +496,7 @@ QgsRasterCalculator::Result QgsRasterCalculator::processCalculationGPU( std::uni
{
cExpression.replace( QStringLiteral( "\"%1\"" ).arg( ref.name ), QStringLiteral( "%1[i]" ).arg( ref.varName ) );
inputArgs.append( QStringLiteral( "__global %1 *%2" )
.arg( ref.typeName )
.arg( ref.varName ) );
.arg( ref.typeName, ref.varName ) );
inputBuffers.push_back( cl::Buffer( ctx, CL_MEM_READ_ONLY, ref.bufferSize, nullptr, nullptr ) );
}

Expand Down Expand Up @@ -490,10 +526,10 @@ QgsRasterCalculator::Result QgsRasterCalculator::processCalculationGPU( std::uni
QStringList inputNoDataCheck;
for ( const auto &ref : inputRefs )
{
inputDesc.append( QStringLiteral( " // %1 = %2" ).arg( ref.varName ).arg( ref.name ) );
inputDesc.append( QStringLiteral( " // %1 = %2" ).arg( ref.varName, ref.name ) );
if ( ref.layer->dataProvider()->sourceHasNoDataValue( ref.band ) )
{
inputNoDataCheck.append( QStringLiteral( "(float) %1[i] == (float) %2" ).arg( ref.varName ).arg( ref.layer->dataProvider()->sourceNoDataValue( ref.band ) ) );
inputNoDataCheck.append( QStringLiteral( "(float) %1[i] == (float) %2" ).arg( ref.varName, QString::number( ref.layer->dataProvider()->sourceNoDataValue( ref.band ), 'g', 10 ) ) );
}
}

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

// qDebug() << programTemplate;
//qDebug() << programTemplate;

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

inputBuffers.clear();

GDALComputeRasterStatistics( outputRasterBand, true, nullptr, nullptr, nullptr, nullptr, GdalProgressCallback, feedback );

}
catch ( cl::Error &e )
{
Expand Down
1 change: 1 addition & 0 deletions src/app/qgisapp.cpp
Expand Up @@ -6807,6 +6807,7 @@ void QgisApp::showRasterCalculator()
QgsProject::instance()->transformContext() );

QProgressDialog p( tr( "Calculating raster expression…" ), tr( "Abort" ), 0, 0 );
p.setWindowTitle( tr( "Raster calculator" ) );
p.setWindowModality( Qt::WindowModal );
p.setMaximum( 100.0 );
QgsFeedback feedback;
Expand Down
48 changes: 41 additions & 7 deletions tests/src/analysis/testqgsrastercalculator.cpp
Expand Up @@ -65,6 +65,8 @@ class TestQgsRasterCalculator : public QObject
void testRasterEntries();
void calcFormulasWithReprojectedLayers();

void testStatistics();

private:

QgsRasterLayer *mpLandsatRasterLayer = nullptr;
Expand Down Expand Up @@ -799,13 +801,13 @@ void TestQgsRasterCalculator::calcFormulasWithReprojectedLayers()
QgsRasterBlock *block = result->dataProvider()->block( 1, extent, 2, 3 );
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 );
qDebug() << "Expected:" << values[0] << values[1] << values[2] << values[3] << values[4] << values[5];
const double epsilon { 0.0001 };
QVERIFY( qgsDoubleNear( block->value( 0, 0 ), static_cast<double>( values[0] ), epsilon ) );
QVERIFY( qgsDoubleNear( block->value( 0, 1 ), static_cast<double>( values[1] ), epsilon ) );
QVERIFY( qgsDoubleNear( block->value( 1, 0 ), static_cast<double>( values[2] ), epsilon ) );
QVERIFY( qgsDoubleNear( block->value( 1, 1 ), static_cast<double>( values[3] ), epsilon ) );
QVERIFY( qgsDoubleNear( block->value( 2, 0 ), static_cast<double>( values[4] ), epsilon ) );
QVERIFY( qgsDoubleNear( block->value( 2, 1 ), static_cast<double>( values[5] ), epsilon ) );
const double epsilon { 0.001 };
QVERIFY2( qgsDoubleNear( block->value( 0, 0 ), static_cast<double>( values[0] ), epsilon ), formula.toUtf8().constData() );
QVERIFY2( qgsDoubleNear( block->value( 0, 1 ), static_cast<double>( values[1] ), epsilon ), formula.toUtf8().constData() );
QVERIFY2( qgsDoubleNear( block->value( 1, 0 ), static_cast<double>( values[2] ), epsilon ), formula.toUtf8().constData() );
QVERIFY2( qgsDoubleNear( block->value( 1, 1 ), static_cast<double>( values[3] ), epsilon ), formula.toUtf8().constData() );
QVERIFY2( qgsDoubleNear( block->value( 2, 0 ), static_cast<double>( values[4] ), epsilon ), formula.toUtf8().constData() );
QVERIFY2( qgsDoubleNear( block->value( 2, 1 ), static_cast<double>( values[5] ), epsilon ), formula.toUtf8().constData() );
delete result;
delete block;
};
Expand Down Expand Up @@ -842,6 +844,38 @@ void TestQgsRasterCalculator::calcFormulasWithReprojectedLayers()

}

void TestQgsRasterCalculator::testStatistics()
{
QgsRasterCalculatorEntry entry1;
entry1.bandNumber = 1;
entry1.raster = mpLandsatRasterLayer;
entry1.ref = QStringLiteral( "landsat@1" );

QTemporaryFile tmpFile;
tmpFile.open(); // fileName is not available until open
QString tmpName = tmpFile.fileName();
tmpFile.close();

QgsCoordinateReferenceSystem crs( QStringLiteral( "EPSG:32633" ) );
QgsRectangle extent( 783235, 3348110, 783350, 3347960 );

QgsRasterCalculator rc( QStringLiteral( "\"landsat@1\" * 2" ),
tmpName,
QStringLiteral( "GTiff" ),
extent, crs, 2, 3, { entry1 },
QgsProject::instance()->transformContext() );
QCOMPARE( static_cast< int >( rc.processCalculation() ), 0 );

//open output file and check stats are there
auto ds = GDALOpenEx( tmpName.toUtf8().constData(), GDAL_OF_RASTER | GDAL_OF_READONLY, nullptr, nullptr, nullptr );
auto band = GDALGetRasterBand( ds, 1 );
double sMin, sMax, sMean, sStdDev;
QCOMPARE( GDALGetRasterStatistics( band, true, false, &sMin, &sMax, &sMean, &sStdDev ), CE_None );
QCOMPARE( sMin, 248.0 );
QCOMPARE( sMax, 250.0 );

}


QGSTEST_MAIN( TestQgsRasterCalculator )
#include "testqgsrastercalculator.moc"

0 comments on commit c1fa7c9

Please sign in to comment.