31 #include <cpl_string.h> 32 #include <gdalwarper.h> 40 : mFormulaString( formulaString )
41 , mOutputFile( outputFile )
42 , mOutputFormat( outputFormat )
43 , mOutputRectangle( outputExtent )
44 , mNumOutputColumns( nOutputColumns )
45 , mNumOutputRows( nOutputRows )
46 , mRasterEntries( rasterEntries )
47 , mTransformContext( transformContext )
55 : mFormulaString( formulaString )
56 , mOutputFile( outputFile )
57 , mOutputFormat( outputFormat )
58 , mOutputRectangle( outputExtent )
59 , mOutputCrs( outputCrs )
60 , mNumOutputColumns( nOutputColumns )
61 , mNumOutputRows( nOutputRows )
62 , mRasterEntries( rasterEntries )
63 , mTransformContext( transformContext )
70 const QgsRectangle &outputExtent,
int nOutputColumns,
int nOutputRows,
const QVector<QgsRasterCalculatorEntry> &rasterEntries )
71 : mFormulaString( formulaString )
72 , mOutputFile( outputFile )
73 , mOutputFormat( outputFormat )
74 , mOutputRectangle( outputExtent )
75 , mNumOutputColumns( nOutputColumns )
76 , mNumOutputRows( nOutputRows )
77 , mRasterEntries( rasterEntries )
80 mOutputCrs = mRasterEntries.at( 0 ).raster->crs();
88 : mFormulaString( formulaString )
89 , mOutputFile( outputFile )
90 , mOutputFormat( outputFormat )
91 , mOutputRectangle( outputExtent )
92 , mOutputCrs( outputCrs )
93 , mNumOutputColumns( nOutputColumns )
94 , mNumOutputRows( nOutputRows )
95 , mRasterEntries( rasterEntries )
113 for (
const auto &entry : qgis::as_const( mRasterEntries ) )
117 mLastError = QObject::tr(
"No raster layer for entry %1" ).arg( entry.ref );
120 if ( entry.bandNumber <= 0 || entry.bandNumber > entry.raster->bandCount() )
122 mLastError = QObject::tr(
"Band number %1 is not valid for entry %2" ).arg( entry.bandNumber ).arg( entry.ref );
130 bool requiresMatrix = ! calcNode->findNodes( QgsRasterCalcNode::Type::tMatrix ).isEmpty();
134 QList<const QgsRasterCalcNode *> nodeList;
137 return processCalculationGPU( std::move( calcNode ), feedback );
142 GDALDriverH outputDriver = openOutputDriver();
145 mLastError = QObject::tr(
"Could not obtain driver for %1" ).arg( mOutputFormat );
150 if ( !outputDataset )
152 mLastError = QObject::tr(
"Could not create output %1" ).arg( mOutputFile );
157 GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
159 float outputNodataValue = -FLT_MAX;
160 GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
164 if ( ! requiresMatrix )
167 std::map<QString, std::unique_ptr<QgsRasterBlock>> inputBlocks;
168 std::map<QString, QgsRasterCalculatorEntry> uniqueRasterEntries;
169 const QList<const QgsRasterCalcNode *> rasterRefNodes = calcNode->findNodes( QgsRasterCalcNode::Type::tRasterRef );
172 QString layerRef( r->toString().remove( 0, 1 ) );
174 if ( ! inputBlocks.count( layerRef ) )
178 if ( ref.ref == layerRef )
180 uniqueRasterEntries[layerRef] = ref;
181 inputBlocks[layerRef ] = qgis::make_unique<QgsRasterBlock>();
188 QMap<QString, QgsRasterBlock * > _rasterData;
190 std::vector<float> castedResult( static_cast<size_t>( mNumOutputColumns ), 0 );
191 auto rowHeight = mOutputRectangle.
height() / mNumOutputRows;
192 for (
size_t row = 0; row < static_cast<size_t>( mNumOutputRows ); ++row )
196 feedback->
setProgress( 100.0 * static_cast< double >( row ) / mNumOutputRows );
210 for (
auto &layerRef : inputBlocks )
219 layerRef.second.reset( proj.
block( ref.
bandNumber, rect, mNumOutputColumns, 1 ) );
228 QgsRasterMatrix resultMatrix( mNumOutputColumns, 1,
nullptr, outputNodataValue );
231 for (
const auto &layerRef : inputBlocks )
233 _rasterData.insert( layerRef.first, layerRef.second.get() );
236 if ( calcNode->calculate( _rasterData, resultMatrix, 0 ) )
238 std::copy( resultMatrix.
data(), resultMatrix.
data() + mNumOutputColumns, castedResult.begin() );
239 if ( GDALRasterIO( outputRasterBand, GF_Write, 0, row, mNumOutputColumns, 1, castedResult.data(), mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
241 QgsDebugMsg( QStringLiteral(
"RasterIO error!" ) );
259 QMap< QString, QgsRasterBlock * > inputBlocks;
260 QVector<QgsRasterCalculatorEntry>::const_iterator it = mRasterEntries.constBegin();
261 for ( ; it != mRasterEntries.constEnd(); ++it )
264 std::unique_ptr< QgsRasterBlock > block;
266 if ( it->raster->crs() != mOutputCrs )
269 proj.
setCrs( it->raster->crs(), mOutputCrs, it->raster->transformContext() );
270 proj.
setInput( it->raster->dataProvider() );
275 block.reset( proj.
block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows, rasterBlockFeedback ) );
278 qDeleteAll( inputBlocks );
284 block.reset( it->raster->dataProvider()->block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows ) );
286 if ( block->isEmpty() )
288 mLastError = QObject::tr(
"Could not allocate required memory for %1" ).arg( it->ref );
289 qDeleteAll( inputBlocks );
292 inputBlocks.insert( it->ref, block.release() );
299 for (
int i = 0; i < mNumOutputRows; ++i )
303 feedback->
setProgress( 100.0 * static_cast< double >( i ) / mNumOutputRows );
311 if ( calcNode->calculate( inputBlocks, resultMatrix, i ) )
313 bool resultIsNumber = resultMatrix.
isNumber();
314 float *calcData =
new float[mNumOutputColumns];
316 for (
int j = 0; j < mNumOutputColumns; ++j )
318 calcData[j] = ( float )( resultIsNumber ? resultMatrix.
number() : resultMatrix.
data()[j] );
322 if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
324 QgsDebugMsg( QStringLiteral(
"RasterIO error!" ) );
331 qDeleteAll( inputBlocks );
346 qDeleteAll( inputBlocks );
364 QString cExpression( calcNode->toString(
true ) );
366 QList<const QgsRasterCalcNode *> nodeList( calcNode->findNodes( QgsRasterCalcNode::Type::tRasterRef ) );
367 QSet<QString> capturedTexts;
368 for (
const auto &r : qgis::as_const( nodeList ) )
370 QString s( r->toString().remove( 0, 1 ) );
372 capturedTexts.insert( s );
389 std::vector<LayerRef> inputRefs;
390 size_t refCounter = 0;
391 for (
const auto &r : capturedTexts )
393 if ( r.startsWith(
'"' ) )
395 QStringList parts( r.split(
'@' ) );
396 if ( parts.count() != 2 )
401 entry.band = parts[1].toInt( &ok );
402 for (
const auto &ref : mRasterEntries )
404 if ( ref.ref == entry.name )
405 entry.layer = ref.raster;
407 if ( !( entry.layer && entry.layer->dataProvider() && ok ) )
409 entry.dataSize = entry.layer->dataProvider()->dataTypeSize( entry.band );
410 switch ( entry.layer->dataProvider()->dataType( entry.band ) )
412 case Qgis::DataType::Byte:
413 entry.typeName = QStringLiteral(
"unsigned char" );
415 case Qgis::DataType::UInt16:
416 entry.typeName = QStringLiteral(
"unsigned int" );
418 case Qgis::DataType::Int16:
419 entry.typeName = QStringLiteral(
"short" );
421 case Qgis::DataType::UInt32:
422 entry.typeName = QStringLiteral(
"unsigned int" );
424 case Qgis::DataType::Int32:
425 entry.typeName = QStringLiteral(
"int" );
427 case Qgis::DataType::Float32:
428 entry.typeName = QStringLiteral(
"float" );
433 case Qgis::DataType::Float64:
434 entry.typeName = QStringLiteral(
"double" );
439 entry.bufferSize = entry.dataSize * mNumOutputColumns;
440 entry.index = refCounter;
441 entry.varName = QStringLiteral(
"input_raster_%1_band_%2" )
444 inputRefs.push_back( entry );
455 std::vector<cl::Buffer> inputBuffers;
456 inputBuffers.reserve( inputRefs.size() );
457 QStringList inputArgs;
458 for (
const auto &ref : inputRefs )
460 cExpression.replace( QStringLiteral(
"\"%1\"" ).arg( ref.name ), QStringLiteral(
"%1[i]" ).arg( ref.varName ) );
461 inputArgs.append( QStringLiteral(
"__global %1 *%2" )
463 .arg( ref.varName ) );
464 inputBuffers.push_back( cl::Buffer( ctx, CL_MEM_READ_ONLY, ref.bufferSize,
nullptr,
nullptr ) );
470 QString programTemplate( R
"CL( 473 // Expression: ##EXPRESSION_ORIGINAL## 474 __kernel void rasterCalculator( ##INPUT## 475 __global float *resultLine 478 // Get the index of the current element 479 const int i = get_global_id(0); 481 resultLine[i] = ##EXPRESSION##; 485 QStringList inputDesc; 486 for (
const auto &ref : inputRefs )
488 inputDesc.append( QStringLiteral(
" // %1 = %2" ).arg( ref.varName ).arg( ref.name ) );
490 programTemplate = programTemplate.replace( QStringLiteral(
"##INPUT_DESC##" ), inputDesc.join(
'\n' ) );
491 programTemplate = programTemplate.replace( QStringLiteral(
"##INPUT##" ), !inputArgs.isEmpty() ? ( inputArgs.join(
',' ).append(
',' ) ) : QChar(
' ' ) );
492 programTemplate = programTemplate.replace( QStringLiteral(
"##EXPRESSION##" ), cExpression );
493 programTemplate = programTemplate.replace( QStringLiteral(
"##EXPRESSION_ORIGINAL##" ), calcNode->toString( ) );
502 Q_ASSERT(
sizeof(
float ) == 4 );
503 std::size_t resultBufferSize( 4 * static_cast<size_t>( mNumOutputColumns ) );
504 cl::Buffer resultLineBuffer( ctx, CL_MEM_WRITE_ONLY,
505 resultBufferSize,
nullptr,
nullptr );
507 auto kernel = cl::Kernel( program,
"rasterCalculator" );
509 for (
unsigned int i = 0; i < inputBuffers.size() ; i++ )
511 kernel.setArg( i, inputBuffers.at( i ) );
513 kernel.setArg( static_cast<unsigned int>( inputBuffers.size() ), resultLineBuffer );
518 GDALDriverH outputDriver = openOutputDriver();
521 mLastError = QObject::tr(
"Could not obtain driver for %1" ).arg( mOutputFormat );
526 if ( !outputDataset )
528 mLastError = QObject::tr(
"Could not create output %1" ).arg( mOutputFile );
535 GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
536 if ( !outputRasterBand )
539 float outputNodataValue = -FLT_MAX;
540 GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
543 std::unique_ptr<QgsRasterBlock> block;
546 auto rowHeight = mOutputRectangle.height() / mNumOutputRows;
547 for (
int line = 0; line < mNumOutputRows; line++ )
556 feedback->
setProgress( 100.0 * static_cast< double >( line ) / mNumOutputRows );
560 for (
const auto &ref : inputRefs )
569 if ( ref.layer->crs() != mOutputCrs )
572 proj.
setCrs( ref.layer->crs(), mOutputCrs, ref.layer->transformContext() );
573 proj.
setInput( ref.layer->dataProvider() );
575 block.reset( proj.
block( ref.band, rect, mNumOutputColumns, 1 ) );
579 block.reset( ref.layer->dataProvider()->block( ref.band, rect, mNumOutputColumns, 1 ) );
586 Q_ASSERT( ref.bufferSize == static_cast<size_t>( block->data().size( ) ) );
587 queue.enqueueWriteBuffer( inputBuffers[ref.index], CL_TRUE, 0,
588 ref.bufferSize, block->bits() );
592 queue.enqueueNDRangeKernel(
595 cl::NDRange( mNumOutputColumns )
599 queue.enqueueReadBuffer( resultLineBuffer, CL_TRUE, 0,
600 resultBufferSize, resultLine.
get() );
605 if ( GDALRasterIO( outputRasterBand, GF_Write, 0, line, mNumOutputColumns, 1, resultLine.
get(), mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
618 inputBuffers.clear();
621 catch ( cl::Error &e )
623 mLastError = e.what();
631 GDALDriverH QgsRasterCalculator::openOutputDriver()
634 GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
652 char **papszOptions =
nullptr;
653 gdal::dataset_unique_ptr outputDataset( GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions ) );
654 if ( !outputDataset )
660 double geotransform[6];
661 outputGeoTransform( geotransform );
662 GDALSetGeoTransform( outputDataset.get(), geotransform );
664 return outputDataset;
667 void QgsRasterCalculator::outputGeoTransform(
double *transform )
const 669 transform[0] = mOutputRectangle.xMinimum();
670 transform[1] = mOutputRectangle.width() / mNumOutputColumns;
672 transform[3] = mOutputRectangle.yMaximum();
674 transform[5] = -mOutputRectangle.height() / mNumOutputRows;
684 QVector<QgsRasterCalculatorEntry> availableEntries;
690 entry.ref = QStringLiteral(
"%1@%2" ).arg( entry.raster->name() ).arg( entry.bandNumber );
694 for (
const auto &ref : qgis::as_const( availableEntries ) )
697 if ( !( entry.raster && ref.raster ) )
700 if ( ref.raster->publicSource() == entry.raster->publicSource() )
702 if ( ref.bandNumber != entry.bandNumber )
712 if ( ref.ref == entry.ref )
715 entry.ref = QStringLiteral(
"%1_%2@%3" ).arg( entry.raster->name() ).arg( i++ ).arg( entry.bandNumber );
723 QMap<QString, QgsMapLayer *>::const_iterator layerIt = layers.constBegin();
724 for ( ; layerIt != layers.constEnd(); ++layerIt )
730 for (
int i = 0; i < rlayer->
bandCount(); ++i )
735 if ( ! uniqueRasterBandIdentifier( entry ) )
737 availableEntries.push_back( entry );
741 return availableEntries;
A rectangle specified with double values.
static QVector< QgsRasterCalculatorEntry > rasterEntries()
Creates a list of raster entries from the current project.
Tiny smart-pointer-like wrapper around CPLMalloc and CPLFree: this is needed because OpenCL C++ API m...
bool isNumber() const
Returns true if matrix is 1x1 (=scalar number)
Result processCalculation(QgsFeedback *feedback=nullptr)
Starts the calculation and writes a new raster.
QString lastError() const
Returns a description of the last error encountered.
QgsRasterLayer * raster
Raster layer associated with entry.
int bandCount() const
Returns the number of bands in this layer.
void cancel()
Tells the internal routines that the current operation should be canceled. This should be run by the ...
Represents a raster layer.
Error allocating memory for result.
void setNodataValue(double d)
void setProgress(double progress)
Sets the current progress for the feedback object.
void canceled()
Internal routines can connect to this signal if they use event loop.
static bool available()
Checks whether a suitable OpenCL platform and device is available on this system and initialize the Q...
static bool supportsRasterCreate(GDALDriverH driver)
Reads whether a driver supports GDALCreate() for raster purposes.
static cl::Context context()
Context factory.
int bandNumber
Band number for entry.
User canceled calculation.
Base class for feedback objects to be used for cancellation of something running in a worker thread...
QgsRasterDataProvider * dataProvider() override
Returns the source data provider.
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...
Error reading input layer.
QString providerType() const
[ data provider interface ] Which provider is being used for this Raster Layer?
void setYMinimum(double y)
Set the minimum y value.
static bool enabled()
Returns true if OpenCL is enabled in the user settings.
Error occurred while performing calculation.
Contains information about the context in which a coordinate transform is executed.
Error creating output data file.
QMap< QString, QgsMapLayer * > mapLayers(const bool validOnly=false) const
Returns a map of all registered layers by layer ID.
QgsRasterBlock * block(int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback=nullptr) override
Read block of data using given extent and size.
Represents an individual raster layer/band number entry within a raster calculation.
Result
Result of the calculation.
void setPrecision(Precision precision)
Full WKT2 string, conforming to ISO 19162:2018 / OGC 18-010, with all possible nodes and new keyword ...
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...
QgsRasterProjector implements approximate projection support for it calculates grid of points in sour...
virtual bool setInput(QgsRasterInterface *input)
Set input.
QgsCoordinateTransformContext transformContext
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.
QgsRasterBlock * block(int bandNo, const QgsRectangle &boundingBox, int width, int height, QgsRasterBlockFeedback *feedback=nullptr) override
Read block of data using given extent and size.
bool isCanceled() const
Tells whether the operation has been canceled already.
void setYMaximum(double y)
Set the maximum y value.
static QgsProject * instance()
Returns the QgsProject singleton instance.
This class represents a coordinate reference system (CRS).
double * data()
Returns data array (but not ownership)
const QgsCoordinateReferenceSystem & outputCrs
static QgsRasterCalcNode * parseRasterCalcString(const QString &str, QString &parserErrorMsg)
double yMaximum() const
Returns the y maximum value (top side of rectangle).
std::unique_ptr< std::remove_pointer< GDALDatasetH >::type, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
Invalid band number for input.
static cl::CommandQueue commandQueue()
Create an OpenCL command queue from the default context.
Feedback object tailored for raster block reading.
QgsCoordinateReferenceSystem crs
QString toWkt(WktVariant variant=WKT1_GDAL, bool multiline=false, int indentationWidth=4) const
Returns a WKT representation of this CRS.
double height() const
Returns the height of the rectangle.
Q_DECL_DEPRECATED void setCrs(const QgsCoordinateReferenceSystem &srcCRS, const QgsCoordinateReferenceSystem &destCRS, int srcDatumTransform=-1, int destDatumTransform=-1)
Sets the source and destination CRS.