31#include <cpl_string.h>
32#include <gdalwarper.h>
44 const char *pszMessage,
47 Q_UNUSED( pszMessage )
49 static double sDfLastComplete = -1.0;
53 if ( sDfLastComplete > dfComplete )
55 if ( sDfLastComplete >= 1.0 )
56 sDfLastComplete = -1.0;
58 sDfLastComplete = dfComplete;
61 if ( std::floor( sDfLastComplete * 10 ) != std::floor( dfComplete * 10 ) )
66 sDfLastComplete = dfComplete;
75 : mFormulaString( formulaString )
76 , mOutputFile( outputFile )
77 , mOutputFormat( outputFormat )
78 , mOutputRectangle( outputExtent )
79 , mNumOutputColumns( nOutputColumns )
80 , mNumOutputRows( nOutputRows )
81 , mRasterEntries( rasterEntries )
82 , mTransformContext( transformContext )
90 : mFormulaString( formulaString )
91 , mOutputFile( outputFile )
92 , mOutputFormat( outputFormat )
93 , mOutputRectangle( outputExtent )
95 , mNumOutputColumns( nOutputColumns )
96 , mNumOutputRows( nOutputRows )
97 , mRasterEntries( rasterEntries )
98 , mTransformContext( transformContext )
105 const QgsRectangle &outputExtent,
int nOutputColumns,
int nOutputRows,
const QVector<QgsRasterCalculatorEntry> &rasterEntries )
106 : mFormulaString( formulaString )
107 , mOutputFile( outputFile )
108 , mOutputFormat( outputFormat )
109 , mOutputRectangle( outputExtent )
110 , mNumOutputColumns( nOutputColumns )
111 , mNumOutputRows( nOutputRows )
112 , mRasterEntries( rasterEntries )
115 mOutputCrs = mRasterEntries.at( 0 ).raster->crs();
123 : mFormulaString( formulaString )
124 , mOutputFile( outputFile )
125 , mOutputFormat( outputFormat )
126 , mOutputRectangle( outputExtent )
128 , mNumOutputColumns( nOutputColumns )
129 , mNumOutputRows( nOutputRows )
130 , mRasterEntries( rasterEntries )
131 , mTransformContext(
QgsProject::instance()->transformContext() )
148 for (
const auto &entry : std::as_const( mRasterEntries ) )
152 mLastError = QObject::tr(
"No raster layer for entry %1" ).arg( entry.ref );
155 if ( entry.bandNumber <= 0 || entry.bandNumber > entry.raster->bandCount() )
157 mLastError = QObject::tr(
"Band number %1 is not valid for entry %2" ).arg( entry.bandNumber ).arg( entry.ref );
171 return processCalculationGPU( std::move( calcNode ), feedback );
176 GDALDriverH outputDriver = openOutputDriver();
179 mLastError = QObject::tr(
"Could not obtain driver for %1" ).arg( mOutputFormat );
184 if ( !outputDataset )
186 mLastError = QObject::tr(
"Could not create output %1" ).arg( mOutputFile );
191 GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
193 float outputNodataValue = -FLT_MAX;
194 GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
198 if ( ! requiresMatrix )
201 std::map<QString, std::unique_ptr<QgsRasterBlock>> inputBlocks;
202 std::map<QString, QgsRasterCalculatorEntry> uniqueRasterEntries;
206 QString layerRef( r->toString().remove( 0, 1 ) );
208 if ( ! inputBlocks.count( layerRef ) )
212 if ( ref.ref == layerRef )
214 uniqueRasterEntries[layerRef] = ref;
215 inputBlocks[layerRef ] = std::make_unique<QgsRasterBlock>();
222 QMap<QString, QgsRasterBlock * > _rasterData;
224 std::vector<float> castedResult(
static_cast<size_t>( mNumOutputColumns ), 0 );
225 auto rowHeight = mOutputRectangle.
height() / mNumOutputRows;
226 for (
size_t row = 0; row < static_cast<size_t>( mNumOutputRows ); ++row )
230 feedback->
setProgress( 100.0 *
static_cast< double >( row ) / mNumOutputRows );
244 for (
auto &layerRef : inputBlocks )
253 layerRef.second.reset( proj.
block( ref.
bandNumber, rect, mNumOutputColumns, 1 ) );
262 QgsRasterMatrix resultMatrix( mNumOutputColumns, 1,
nullptr, outputNodataValue );
265 for (
const auto &layerRef : inputBlocks )
267 _rasterData.insert( layerRef.first, layerRef.second.get() );
270 if ( calcNode->calculate( _rasterData, resultMatrix, 0 ) )
272 std::copy( resultMatrix.
data(), resultMatrix.
data() + mNumOutputColumns, castedResult.begin() );
273 if ( GDALRasterIO( outputRasterBand, GF_Write, 0, row, mNumOutputColumns, 1, castedResult.data(), mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
293 QMap< QString, QgsRasterBlock * > inputBlocks;
294 QVector<QgsRasterCalculatorEntry>::const_iterator it = mRasterEntries.constBegin();
295 for ( ; it != mRasterEntries.constEnd(); ++it )
298 std::unique_ptr< QgsRasterBlock > block;
300 if ( it->raster->crs() != mOutputCrs )
303 proj.
setCrs( it->raster->crs(), mOutputCrs, it->raster->transformContext() );
304 proj.
setInput( it->raster->dataProvider() );
309 block.reset( proj.
block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows, rasterBlockFeedback ) );
312 qDeleteAll( inputBlocks );
318 block.reset( it->raster->dataProvider()->block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows ) );
320 if ( block->isEmpty() )
322 mLastError = QObject::tr(
"Could not allocate required memory for %1" ).arg( it->ref );
323 qDeleteAll( inputBlocks );
326 inputBlocks.insert( it->ref, block.release() );
333 for (
int i = 0; i < mNumOutputRows; ++i )
337 feedback->
setProgress( 100.0 *
static_cast< double >( i ) / mNumOutputRows );
345 if ( calcNode->calculate( inputBlocks, resultMatrix, i ) )
347 bool resultIsNumber = resultMatrix.
isNumber();
348 float *calcData =
new float[mNumOutputColumns];
350 for (
int j = 0; j < mNumOutputColumns; ++j )
352 calcData[j] = ( float )( resultIsNumber ? resultMatrix.
number() : resultMatrix.
data()[j] );
356 if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
365 qDeleteAll( inputBlocks );
380 qDeleteAll( inputBlocks );
392 GDALComputeRasterStatistics( outputRasterBand,
true,
nullptr,
nullptr,
nullptr,
nullptr,
GdalProgressCallback, feedback );
401 QString cExpression( calcNode->toString(
true ) );
404 QSet<QString> capturedTexts;
405 for (
const auto &r : std::as_const( nodeList ) )
407 QString s( r->toString().remove( 0, 1 ) );
409 capturedTexts.insert( s );
426 std::vector<LayerRef> inputRefs;
427 size_t refCounter = 0;
428 for (
const auto &r : capturedTexts )
430 if ( r.startsWith(
'"' ) )
432 QStringList parts( r.split(
'@' ) );
433 if ( parts.count() != 2 )
438 entry.band = parts[1].toInt( &ok );
439 for (
const auto &ref : std::as_const( mRasterEntries ) )
441 if ( ref.ref == entry.name )
442 entry.layer = ref.raster;
444 if ( !( entry.layer && entry.layer->dataProvider() && ok ) )
446 entry.dataSize = entry.layer->dataProvider()->dataTypeSize( entry.band );
447 switch ( entry.layer->dataProvider()->dataType( entry.band ) )
450 entry.typeName = QStringLiteral(
"unsigned char" );
453 entry.typeName = QStringLiteral(
"signed char" );
456 entry.typeName = QStringLiteral(
"unsigned int" );
459 entry.typeName = QStringLiteral(
"short" );
462 entry.typeName = QStringLiteral(
"unsigned int" );
465 entry.typeName = QStringLiteral(
"int" );
468 entry.typeName = QStringLiteral(
"float" );
474 entry.typeName = QStringLiteral(
"double" );
485 entry.bufferSize = entry.dataSize * mNumOutputColumns;
486 entry.index = refCounter;
487 entry.varName = QStringLiteral(
"input_raster_%1_band_%2" )
490 inputRefs.push_back( entry );
501 std::vector<cl::Buffer> inputBuffers;
502 inputBuffers.reserve( inputRefs.size() );
503 QStringList inputArgs;
504 for (
const auto &ref : inputRefs )
506 cExpression.replace( QStringLiteral(
"\"%1\"" ).arg( ref.name ), QStringLiteral(
"%1[i]" ).arg( ref.varName ) );
507 inputArgs.append( QStringLiteral(
"__global %1 *%2" )
508 .arg( ref.typeName, ref.varName ) );
509 inputBuffers.push_back( cl::Buffer( ctx, CL_MEM_READ_ONLY, ref.bufferSize,
nullptr,
nullptr ) );
515 QString programTemplate( R
"CL(
518 // Expression: ##EXPRESSION_ORIGINAL##
519 __kernel void rasterCalculator( ##INPUT##
520 __global float *resultLine
523 // Get the index of the current element
524 const int i = get_global_id(0);
525 // Check for nodata in input
526 if ( ##INPUT_NODATA_CHECK## )
527 resultLine[i] = -FLT_MAX;
530 resultLine[i] = ##EXPRESSION##;
534 QStringList inputDesc;
535 QStringList inputNoDataCheck;
536 for (
const auto &ref : inputRefs )
538 inputDesc.append( QStringLiteral(
" // %1 = %2" ).arg( ref.varName, ref.name ) );
539 if ( ref.layer->dataProvider()->sourceHasNoDataValue( ref.band ) )
541 inputNoDataCheck.append( QStringLiteral(
"(float) %1[i] == (float) %2" ).arg( ref.varName, QString::number( ref.layer->dataProvider()->sourceNoDataValue( ref.band ),
'g', 10 ) ) );
545 programTemplate = programTemplate.replace( QLatin1String(
"##INPUT_NODATA_CHECK##" ), inputNoDataCheck.isEmpty() ? QStringLiteral(
"false" ) : inputNoDataCheck.join( QLatin1String(
" || " ) ) );
546 programTemplate = programTemplate.replace( QLatin1String(
"##INPUT_DESC##" ), inputDesc.join(
'\n' ) );
547 programTemplate = programTemplate.replace( QLatin1String(
"##INPUT##" ), !inputArgs.isEmpty() ? ( inputArgs.join(
',' ).append(
',' ) ) : QChar(
' ' ) );
548 programTemplate = programTemplate.replace( QLatin1String(
"##EXPRESSION##" ), cExpression );
549 programTemplate = programTemplate.replace( QLatin1String(
"##EXPRESSION_ORIGINAL##" ), calcNode->toString( ) );
558 Q_ASSERT(
sizeof(
float ) == 4 );
559 std::size_t resultBufferSize( 4 *
static_cast<size_t>( mNumOutputColumns ) );
560 cl::Buffer resultLineBuffer( ctx, CL_MEM_WRITE_ONLY,
561 resultBufferSize,
nullptr,
nullptr );
563 auto kernel = cl::Kernel( program,
"rasterCalculator" );
565 for (
unsigned int i = 0; i < inputBuffers.size() ; i++ )
567 kernel.setArg( i, inputBuffers.at( i ) );
569 kernel.setArg(
static_cast<unsigned int>( inputBuffers.size() ), resultLineBuffer );
574 GDALDriverH outputDriver = openOutputDriver();
577 mLastError = QObject::tr(
"Could not obtain driver for %1" ).arg( mOutputFormat );
582 if ( !outputDataset )
584 mLastError = QObject::tr(
"Could not create output %1" ).arg( mOutputFile );
591 GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
592 if ( !outputRasterBand )
595 const float outputNodataValue = -FLT_MAX;
596 GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
599 std::unique_ptr<QgsRasterBlock> block;
602 auto rowHeight = mOutputRectangle.
height() / mNumOutputRows;
603 for (
int line = 0; line < mNumOutputRows; line++ )
612 feedback->
setProgress( 100.0 *
static_cast< double >( line ) / mNumOutputRows );
616 for (
const auto &ref : inputRefs )
620 rect.setYMaximum( rect.yMaximum() - rowHeight * line );
621 rect.setYMinimum( rect.yMaximum() - rowHeight );
625 if ( ref.layer->crs() != mOutputCrs )
628 proj.
setCrs( ref.layer->crs(), mOutputCrs, ref.layer->transformContext() );
629 proj.
setInput( ref.layer->dataProvider() );
631 block.reset( proj.
block( ref.band, rect, mNumOutputColumns, 1 ) );
635 block.reset( ref.layer->dataProvider()->block( ref.band, rect, mNumOutputColumns, 1 ) );
642 Q_ASSERT( ref.bufferSize ==
static_cast<size_t>( block->data().size( ) ) );
643 queue.enqueueWriteBuffer( inputBuffers[ref.index], CL_TRUE, 0,
644 ref.bufferSize, block->bits() );
648 queue.enqueueNDRangeKernel(
651 cl::NDRange( mNumOutputColumns )
655 queue.enqueueReadBuffer( resultLineBuffer, CL_TRUE, 0,
656 resultBufferSize, resultLine.get() );
661 if ( GDALRasterIO( outputRasterBand, GF_Write, 0, line, mNumOutputColumns, 1, resultLine.get(), mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
674 inputBuffers.clear();
676 GDALComputeRasterStatistics( outputRasterBand,
true,
nullptr,
nullptr,
nullptr,
nullptr,
GdalProgressCallback, feedback );
679 catch ( cl::Error &e )
681 mLastError = e.what();
689GDALDriverH QgsRasterCalculator::openOutputDriver()
692 GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
710 char **papszOptions =
nullptr;
711 gdal::dataset_unique_ptr outputDataset( GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions ) );
712 if ( !outputDataset )
718 double geotransform[6];
719 outputGeoTransform( geotransform );
720 GDALSetGeoTransform( outputDataset.get(), geotransform );
722 return outputDataset;
725void QgsRasterCalculator::outputGeoTransform(
double *transform )
const
727 transform[0] = mOutputRectangle.
xMinimum();
728 transform[1] = mOutputRectangle.
width() / mNumOutputColumns;
730 transform[3] = mOutputRectangle.
yMaximum();
732 transform[5] = -mOutputRectangle.
height() / mNumOutputRows;
742 QVector<QgsRasterCalculatorEntry> availableEntries;
748 entry.ref = QStringLiteral(
"%1@%2" ).arg( entry.raster->name() ).arg( entry.bandNumber );
752 for (
const auto &
ref : std::as_const( availableEntries ) )
755 if ( !( entry.raster &&
ref.raster ) )
758 if (
ref.raster->publicSource() == entry.raster->publicSource() )
760 if (
ref.bandNumber != entry.bandNumber )
770 if (
ref.ref == entry.ref )
773 entry.ref = QStringLiteral(
"%1_%2@%3" ).arg( entry.raster->name() ).arg( i++ ).arg( entry.bandNumber );
781 QMap<QString, QgsMapLayer *>::const_iterator layerIt = layers.constBegin();
782 for ( ; layerIt != layers.constEnd(); ++layerIt )
784 QgsRasterLayer *rlayer = qobject_cast<QgsRasterLayer *>( layerIt.value() );
788 for (
int i = 0; i < rlayer->
bandCount(); ++i )
793 if ( ! uniqueRasterBandIdentifier( entry ) )
795 availableEntries.push_back( entry );
799 return availableEntries;
@ Size
Original data source size (and thus resolution) is known, it is not always available,...
@ Float32
Thirty two bit floating point (float)
@ CFloat64
Complex Float64.
@ Int16
Sixteen bit signed integer (qint16)
@ ARGB32_Premultiplied
Color, alpha, red, green, blue, 4 bytes the same as QImage::Format_ARGB32_Premultiplied.
@ Int8
Eight bit signed integer (qint8) (added in QGIS 3.30)
@ UInt16
Sixteen bit unsigned integer (quint16)
@ Byte
Eight bit unsigned integer (quint8)
@ UnknownDataType
Unknown or unspecified type.
@ ARGB32
Color, alpha, red, green, blue, 4 bytes the same as QImage::Format_ARGB32.
@ Int32
Thirty two bit signed integer (qint32)
@ Float64
Sixty four bit floating point (double)
@ CFloat32
Complex Float32.
@ UInt32
Thirty two bit unsigned integer (quint32)
@ PreferredGdal
Preferred format for conversion of CRS to WKT for use with the GDAL library.
This class represents a coordinate reference system (CRS).
QString toWkt(Qgis::CrsWktVariant variant=Qgis::CrsWktVariant::Wkt1Gdal, bool multiline=false, int indentationWidth=4) const
Returns a WKT representation of this CRS.
Contains information about the context in which a coordinate transform is executed.
Base class for feedback objects to be used for cancellation of something running in a worker thread.
bool isCanceled() const
Tells whether the operation has been canceled already.
void canceled()
Internal routines can connect to this signal if they use event loop.
void cancel()
Tells the internal routines that the current operation should be canceled. This should be run by the ...
void setProgress(double progress)
Sets the current progress for the feedback object.
static bool supportsRasterCreate(GDALDriverH driver)
Reads whether a driver supports GDALCreate() for raster purposes.
QgsCoordinateReferenceSystem crs
static Q_DECL_DEPRECATED cl::Program buildProgram(const cl::Context &context, const QString &source, ExceptionBehavior exceptionBehavior=Catch)
Build the program from source in the given context and depending on exceptionBehavior can throw or ca...
static cl::Context context()
Context factory.
static bool enabled()
Returns true if OpenCL is enabled in the user settings.
static bool available()
Checks whether a suitable OpenCL platform and device is available on this system and initialize the Q...
static cl::CommandQueue commandQueue()
Create an OpenCL command queue from the default context.
@ Throw
Write errors in the message log and re-throw exceptions.
Encapsulates a QGIS project, including sets of map layers and their styles, layouts,...
static QgsProject * instance()
Returns the QgsProject singleton instance.
QgsCoordinateTransformContext transformContext
QMap< QString, QgsMapLayer * > mapLayers(const bool validOnly=false) const
Returns a map of all registered layers by layer ID.
Feedback object tailored for raster block reading.
Represents a node in a raster calculator.
static QgsRasterCalcNode * parseRasterCalcString(const QString &str, QString &parserErrorMsg)
Represents an individual raster layer/band number entry within a raster calculation.
QgsRasterLayer * raster
Raster layer associated with entry.
static QVector< QgsRasterCalculatorEntry > rasterEntries()
Creates a list of raster entries from the current project.
int bandNumber
Band number for entry.
QString ref
Name of entry.
QgsRasterCalculator(const QString &formulaString, const QString &outputFile, const QString &outputFormat, const QgsRectangle &outputExtent, int nOutputColumns, int nOutputRows, const QVector< QgsRasterCalculatorEntry > &rasterEntries, const QgsCoordinateTransformContext &transformContext)
QgsRasterCalculator constructor.
QString lastError() const
Returns a description of the last error encountered.
Result processCalculation(QgsFeedback *feedback=nullptr)
Starts the calculation and writes a new raster.
Result
Result of the calculation.
@ InputLayerError
Error reading input layer.
@ ParserError
Error parsing formula.
@ CalculationError
Error occurred while performing calculation.
@ MemoryError
Error allocating memory for result.
@ Canceled
User canceled calculation.
@ Success
Calculation successful.
@ BandError
Invalid band number for input.
@ CreateOutputError
Error creating output data file.
QgsRasterBlock * block(int bandNo, const QgsRectangle &boundingBox, int width, int height, QgsRasterBlockFeedback *feedback=nullptr) override
Read block of data using given extent and size.
virtual Qgis::RasterInterfaceCapabilities capabilities() const
Returns the capabilities supported by the interface.
virtual bool setInput(QgsRasterInterface *input)
Set input.
Represents a raster layer.
int bandCount() const
Returns the number of bands in this layer.
QgsRasterDataProvider * dataProvider() override
Returns the source data provider.
Represents a matrix in a raster calculator operation.
bool isNumber() const
Returns true if matrix is 1x1 (=scalar number)
void setNodataValue(double d)
double * data()
Returns data array (but not ownership)
Implements approximate projection support for optimised raster transformation.
void setPrecision(Precision precision)
QgsRasterBlock * block(int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback=nullptr) override
Read block of data using given extent and size.
@ Exact
Exact, precise but slow.
Q_DECL_DEPRECATED void setCrs(const QgsCoordinateReferenceSystem &srcCRS, const QgsCoordinateReferenceSystem &destCRS, int srcDatumTransform=-1, int destDatumTransform=-1)
Sets the source and destination CRS.
A rectangle specified with double values.
double xMinimum() const
Returns the x minimum value (left side of rectangle).
void setYMinimum(double y)
Set the minimum y value.
double width() const
Returns the width of the rectangle.
double yMaximum() const
Returns the y maximum value (top side of rectangle).
void setYMaximum(double y)
Set the maximum y value.
double height() const
Returns the height of the rectangle.
void CORE_EXPORT fast_delete_and_close(dataset_unique_ptr &dataset, GDALDriverH driver, const QString &path)
Performs a fast close of an unwanted GDAL dataset handle by deleting the underlying data store.
std::unique_ptr< std::remove_pointer< GDALDatasetH >::type, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
#define QgsDebugError(str)
int CPL_STDCALL GdalProgressCallback(double dfComplete, const char *pszMessage, void *pProgressArg)
const QgsCoordinateReferenceSystem & outputCrs
Tiny smart-pointer-like wrapper around CPLMalloc and CPLFree: this is needed because OpenCL C++ API m...