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 );
129 QList<const QgsRasterCalcNode *> nodeList;
131 return processCalculationGPU( std::move( calcNode ), feedback );
135 GDALDriverH outputDriver = openOutputDriver();
138 mLastError = QObject::tr(
"Could not obtain driver for %1" ).arg( mOutputFormat );
143 if ( !outputDataset )
145 mLastError = QObject::tr(
"Could not create output %1" ).arg( mOutputFile );
149 GDALSetProjection( outputDataset.get(), mOutputCrs.
toWkt().toLocal8Bit().data() );
150 GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
152 float outputNodataValue = -FLT_MAX;
153 GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
158 bool requiresMatrix = ! calcNode->findNodes( QgsRasterCalcNode::Type::tMatrix ).isEmpty();
161 if ( ! requiresMatrix )
164 std::map<QString, std::unique_ptr<QgsRasterBlock>> inputBlocks;
165 std::map<QString, QgsRasterCalculatorEntry> uniqueRasterEntries;
166 for (
const auto &r : calcNode->findNodes( QgsRasterCalcNode::Type::tRasterRef ) )
168 QString layerRef( r->toString().remove( 0, 1 ) );
170 if ( ! inputBlocks.count( layerRef ) )
172 for (
const auto &ref : mRasterEntries )
174 if ( ref.ref == layerRef )
176 uniqueRasterEntries[layerRef] = ref;
177 inputBlocks[layerRef ] = qgis::make_unique<QgsRasterBlock>();
184 QMap<QString, QgsRasterBlock * > _rasterData;
186 std::vector<float> castedResult( static_cast<size_t>( mNumOutputColumns ), 0 );
187 auto rowHeight = mOutputRectangle.
height() / mNumOutputRows;
188 for (
size_t row = 0; row < static_cast<size_t>( mNumOutputRows ); ++row )
192 feedback->
setProgress( 100.0 * static_cast< double >( row ) / mNumOutputRows );
206 for (
auto &layerRef : inputBlocks )
209 if ( uniqueRasterEntries[layerRef.first].raster->crs() != mOutputCrs )
215 layerRef.second.reset( proj.
block( ref.
bandNumber, rect, mNumOutputColumns, 1 ) );
227 for (
const auto &layerRef : inputBlocks )
229 _rasterData.insert( layerRef.first, inputBlocks[layerRef.first].get() );
232 if ( calcNode->calculate( _rasterData, resultMatrix, 0 ) )
234 std::copy( resultMatrix.
data(), resultMatrix.
data() + mNumOutputColumns, castedResult.begin() );
235 if ( GDALRasterIO( outputRasterBand, GF_Write, 0, row, mNumOutputColumns, 1, castedResult.data(), mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
237 QgsDebugMsg( QStringLiteral(
"RasterIO error!" ) );
249 QMap< QString, QgsRasterBlock * > inputBlocks;
250 QVector<QgsRasterCalculatorEntry>::const_iterator it = mRasterEntries.constBegin();
251 for ( ; it != mRasterEntries.constEnd(); ++it )
254 std::unique_ptr< QgsRasterBlock > block;
256 if ( it->raster->crs() != mOutputCrs )
259 proj.
setCrs( it->raster->crs(), mOutputCrs, it->raster->transformContext() );
260 proj.
setInput( it->raster->dataProvider() );
265 block.reset( proj.
block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows, rasterBlockFeedback ) );
268 qDeleteAll( inputBlocks );
274 block.reset( it->raster->dataProvider()->block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows ) );
276 if ( block->isEmpty() )
278 mLastError = QObject::tr(
"Could not allocate required memory for %1" ).arg( it->ref );
279 qDeleteAll( inputBlocks );
282 inputBlocks.insert( it->ref, block.release() );
289 for (
int i = 0; i < mNumOutputRows; ++i )
293 feedback->
setProgress( 100.0 * static_cast< double >( i ) / mNumOutputRows );
301 if ( calcNode->calculate( inputBlocks, resultMatrix, i ) )
303 bool resultIsNumber = resultMatrix.
isNumber();
304 float *calcData =
new float[mNumOutputColumns];
306 for (
int j = 0; j < mNumOutputColumns; ++j )
308 calcData[j] = ( float )( resultIsNumber ? resultMatrix.
number() : resultMatrix.
data()[j] );
312 if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
314 QgsDebugMsg( QStringLiteral(
"RasterIO error!" ) );
329 qDeleteAll( inputBlocks );
347 QString cExpression( calcNode->toString(
true ) );
349 QList<const QgsRasterCalcNode *> nodeList( calcNode->findNodes( QgsRasterCalcNode::Type::tRasterRef ) );
350 QSet<QString> capturedTexts;
351 for (
const auto &r : qgis::as_const( nodeList ) )
353 QString s( r->toString().remove( 0, 1 ) );
355 capturedTexts.insert( s );
372 std::vector<LayerRef> inputRefs;
373 size_t refCounter = 0;
374 for (
const auto &r : capturedTexts )
376 if ( r.startsWith(
'"' ) )
378 QStringList parts( r.split(
'@' ) );
379 if ( parts.count() != 2 )
384 entry.band = parts[1].toInt( &ok );
385 for (
const auto &ref : mRasterEntries )
387 if ( ref.ref == entry.name )
388 entry.layer = ref.raster;
390 if ( !( entry.layer && entry.layer->dataProvider() && ok ) )
392 entry.dataSize = entry.layer->dataProvider()->dataTypeSize( entry.band );
393 switch ( entry.layer->dataProvider()->dataType( entry.band ) )
395 case Qgis::DataType::Byte:
396 entry.typeName = QStringLiteral(
"unsigned char" );
398 case Qgis::DataType::UInt16:
399 entry.typeName = QStringLiteral(
"unsigned int" );
401 case Qgis::DataType::Int16:
402 entry.typeName = QStringLiteral(
"short" );
404 case Qgis::DataType::UInt32:
405 entry.typeName = QStringLiteral(
"unsigned int" );
407 case Qgis::DataType::Int32:
408 entry.typeName = QStringLiteral(
"int" );
410 case Qgis::DataType::Float32:
411 entry.typeName = QStringLiteral(
"float" );
416 case Qgis::DataType::Float64:
417 entry.typeName = QStringLiteral(
"double" );
422 entry.bufferSize = entry.dataSize * mNumOutputColumns;
423 entry.index = refCounter;
424 entry.varName = QStringLiteral(
"input_raster_%1_band_%2" )
427 inputRefs.push_back( entry );
435 std::vector<cl::Buffer> inputBuffers;
436 inputBuffers.reserve( inputRefs.size() );
437 QStringList inputArgs;
438 for (
const auto &ref : inputRefs )
440 cExpression.replace( QStringLiteral(
"\"%1\"" ).arg( ref.name ), QStringLiteral(
"%1[i]" ).arg( ref.varName ) );
441 inputArgs.append( QStringLiteral(
"__global %1 *%2" )
443 .arg( ref.varName ) );
444 inputBuffers.push_back( cl::Buffer( ctx, CL_MEM_READ_ONLY, ref.bufferSize,
nullptr,
nullptr ) );
450 QString programTemplate( R
"CL( 453 // Expression: ##EXPRESSION_ORIGINAL## 454 __kernel void rasterCalculator( ##INPUT## 455 __global float *resultLine 458 // Get the index of the current element 459 const int i = get_global_id(0); 461 resultLine[i] = ##EXPRESSION##; 465 QStringList inputDesc; 466 for (
const auto &ref : inputRefs )
468 inputDesc.append( QStringLiteral(
" // %1 = %2" ).arg( ref.varName ).arg( ref.name ) );
470 programTemplate = programTemplate.replace( QStringLiteral(
"##INPUT_DESC##" ), inputDesc.join(
'\n' ) );
471 programTemplate = programTemplate.replace( QStringLiteral(
"##INPUT##" ), !inputArgs.isEmpty() ? ( inputArgs.join(
',' ).append(
',' ) ) : QChar(
' ' ) );
472 programTemplate = programTemplate.replace( QStringLiteral(
"##EXPRESSION##" ), cExpression );
473 programTemplate = programTemplate.replace( QStringLiteral(
"##EXPRESSION_ORIGINAL##" ), calcNode->toString( ) );
482 Q_ASSERT(
sizeof(
float ) == 4 );
483 std::size_t resultBufferSize( 4 * static_cast<size_t>( mNumOutputColumns ) );
484 cl::Buffer resultLineBuffer( ctx, CL_MEM_WRITE_ONLY,
485 resultBufferSize,
nullptr,
nullptr );
487 auto kernel = cl::Kernel( program,
"rasterCalculator" );
489 for (
unsigned int i = 0; i < inputBuffers.size() ; i++ )
491 kernel.setArg( i, inputBuffers.at( i ) );
493 kernel.setArg( static_cast<unsigned int>( inputBuffers.size() ), resultLineBuffer );
498 GDALDriverH outputDriver = openOutputDriver();
501 mLastError = QObject::tr(
"Could not obtain driver for %1" ).arg( mOutputFormat );
506 if ( !outputDataset )
508 mLastError = QObject::tr(
"Could not create output %1" ).arg( mOutputFile );
512 GDALSetProjection( outputDataset.get(), mOutputCrs.toWkt().toLocal8Bit().data() );
514 GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
515 if ( !outputRasterBand )
519 std::unique_ptr<QgsRasterBlock> block;
522 auto rowHeight = mOutputRectangle.height() / mNumOutputRows;
523 for (
int line = 0; line < mNumOutputRows; line++ )
532 feedback->
setProgress( 100.0 * static_cast< double >( line ) / mNumOutputRows );
536 for (
const auto &ref : inputRefs )
545 if ( ref.layer->crs() != mOutputCrs )
548 proj.
setCrs( ref.layer->crs(), mOutputCrs, ref.layer->transformContext() );
549 proj.
setInput( ref.layer->dataProvider() );
551 block.reset( proj.
block( ref.band, rect, mNumOutputColumns, 1 ) );
555 block.reset( ref.layer->dataProvider()->block( ref.band, rect, mNumOutputColumns, 1 ) );
562 Q_ASSERT( ref.bufferSize == static_cast<size_t>( block->data().size( ) ) );
563 queue.enqueueWriteBuffer( inputBuffers[ref.index], CL_TRUE, 0,
564 ref.bufferSize, block->bits() );
568 queue.enqueueNDRangeKernel(
571 cl::NDRange( mNumOutputColumns )
575 queue.enqueueReadBuffer( resultLineBuffer, CL_TRUE, 0,
576 resultBufferSize, resultLine.
get() );
581 if ( GDALRasterIO( outputRasterBand, GF_Write, 0, line, mNumOutputColumns, 1, resultLine.
get(), mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
594 inputBuffers.clear();
600 GDALDriverH QgsRasterCalculator::openOutputDriver()
603 GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
621 char **papszOptions =
nullptr;
622 gdal::dataset_unique_ptr outputDataset( GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions ) );
623 if ( !outputDataset )
629 double geotransform[6];
630 outputGeoTransform( geotransform );
631 GDALSetGeoTransform( outputDataset.get(), geotransform );
633 return outputDataset;
636 void QgsRasterCalculator::outputGeoTransform(
double *transform )
const 638 transform[0] = mOutputRectangle.xMinimum();
639 transform[1] = mOutputRectangle.width() / mNumOutputColumns;
641 transform[3] = mOutputRectangle.yMaximum();
643 transform[5] = -mOutputRectangle.height() / mNumOutputRows;
653 QVector<QgsRasterCalculatorEntry> availableEntries;
659 entry.ref = QStringLiteral(
"%1@%2" ).arg( entry.raster->name() ).arg( entry.bandNumber );
663 for (
const auto &ref : qgis::as_const( availableEntries ) )
666 if ( !( entry.raster && ref.raster ) )
669 if ( ref.raster->publicSource() == entry.raster->publicSource() )
671 if ( ref.bandNumber != entry.bandNumber )
681 if ( ref.ref == entry.ref )
684 entry.ref = QStringLiteral(
"%1_%2@%3" ).arg( entry.raster->name() ).arg( i++ ).arg( entry.bandNumber );
692 QMap<QString, QgsMapLayer *>::const_iterator layerIt = layers.constBegin();
693 for ( ; layerIt != layers.constEnd(); ++layerIt )
699 for (
int i = 0; i < rlayer->
bandCount(); ++i )
704 if ( ! uniqueRasterBandIdentifier( entry ) )
706 availableEntries.push_back( entry );
710 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 ...
This class provides qgis with the ability to render raster datasets onto the mapcanvas.
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.
virtual QString name() const =0
Returns a provider name.
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 layer's data provider, it may be nullptr.
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.
void setYMinimum(double y)
Set the minimum y value.
static bool enabled()
Returns true if OpenCL is enabled in the user settings.
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)
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).
QString toWkt() const
Returns a WKT representation of this 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
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.