QGIS API Documentation  3.27.0-Master (e113457133)
qgsrastercalculator.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsrastercalculator.cpp - description
3  -----------------------
4  begin : September 28th, 2010
5  copyright : (C) 2010 by Marco Hugentobler
6  email : marco dot hugentobler at sourcepole dot ch
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
18 #include "qgsgdalutils.h"
19 #include "qgsrastercalculator.h"
20 #include "qgsrasterdataprovider.h"
21 #include "qgsrasterinterface.h"
22 #include "qgsrasterlayer.h"
23 #include "qgsrastermatrix.h"
24 #include "qgsrasterprojector.h"
25 #include "qgsfeedback.h"
26 #include "qgsogrutils.h"
27 #include "qgsproject.h"
28 
29 #include <QFile>
30 
31 #include <cpl_string.h>
32 #include <gdalwarper.h>
33 
34 #ifdef HAVE_OPENCL
35 #include "qgsopenclutils.h"
36 #include "qgsgdalutils.h"
37 #endif
38 
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 
74 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 )
75  : mFormulaString( formulaString )
76  , mOutputFile( outputFile )
77  , mOutputFormat( outputFormat )
78  , mOutputRectangle( outputExtent )
79  , mNumOutputColumns( nOutputColumns )
80  , mNumOutputRows( nOutputRows )
81  , mRasterEntries( rasterEntries )
82  , mTransformContext( transformContext )
83 {
84 
85 }
86 
87 QgsRasterCalculator::QgsRasterCalculator( const QString &formulaString, const QString &outputFile, const QString &outputFormat,
88  const QgsRectangle &outputExtent, const QgsCoordinateReferenceSystem &outputCrs, int nOutputColumns, int nOutputRows,
89  const QVector<QgsRasterCalculatorEntry> &rasterEntries, const QgsCoordinateTransformContext &transformContext )
90  : mFormulaString( formulaString )
91  , mOutputFile( outputFile )
92  , mOutputFormat( outputFormat )
93  , mOutputRectangle( outputExtent )
94  , mOutputCrs( outputCrs )
95  , mNumOutputColumns( nOutputColumns )
96  , mNumOutputRows( nOutputRows )
97  , mRasterEntries( rasterEntries )
98  , mTransformContext( transformContext )
99 {
100 
101 }
102 
103 // Deprecated!
104 QgsRasterCalculator::QgsRasterCalculator( const QString &formulaString, const QString &outputFile, const QString &outputFormat,
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 )
113 {
114  //default to first layer's crs
115  mOutputCrs = mRasterEntries.at( 0 ).raster->crs();
116  mTransformContext = QgsProject::instance()->transformContext();
117 }
118 
119 
120 // Deprecated!
121 QgsRasterCalculator::QgsRasterCalculator( const QString &formulaString, const QString &outputFile, const QString &outputFormat,
122  const QgsRectangle &outputExtent, const QgsCoordinateReferenceSystem &outputCrs, int nOutputColumns, int nOutputRows, const QVector<QgsRasterCalculatorEntry> &rasterEntries )
123  : mFormulaString( formulaString )
124  , mOutputFile( outputFile )
125  , mOutputFormat( outputFormat )
126  , mOutputRectangle( outputExtent )
127  , mOutputCrs( outputCrs )
128  , mNumOutputColumns( nOutputColumns )
129  , mNumOutputRows( nOutputRows )
130  , mRasterEntries( rasterEntries )
131  , mTransformContext( QgsProject::instance()->transformContext() )
132 {
133 }
134 
136 {
137  mLastError.clear();
138 
139  //prepare search string / tree
140  std::unique_ptr< QgsRasterCalcNode > calcNode( QgsRasterCalcNode::parseRasterCalcString( mFormulaString, mLastError ) );
141  if ( !calcNode )
142  {
143  //error
144  return ParserError;
145  }
146 
147  // Check input layers and bands
148  for ( const auto &entry : std::as_const( mRasterEntries ) )
149  {
150  if ( !entry.raster ) // no raster layer in entry
151  {
152  mLastError = QObject::tr( "No raster layer for entry %1" ).arg( entry.ref );
153  return InputLayerError;
154  }
155  if ( entry.bandNumber <= 0 || entry.bandNumber > entry.raster->bandCount() )
156  {
157  mLastError = QObject::tr( "Band number %1 is not valid for entry %2" ).arg( entry.bandNumber ).arg( entry.ref );
158  return BandError;
159  }
160  }
161 
162  // Check if we need to read the raster as a whole (which is memory inefficient
163  // and not interruptible by the user) by checking if any raster matrix nodes are
164  // in the expression
165  bool requiresMatrix = ! calcNode->findNodes( QgsRasterCalcNode::Type::tMatrix ).isEmpty();
166 
167 #ifdef HAVE_OPENCL
168  // Check for matrix nodes, GPU implementation does not support them
169  if ( QgsOpenClUtils::enabled() && QgsOpenClUtils::available() && ! requiresMatrix )
170  {
171  return processCalculationGPU( std::move( calcNode ), feedback );
172  }
173 #endif
174 
175  //open output dataset for writing
176  GDALDriverH outputDriver = openOutputDriver();
177  if ( !outputDriver )
178  {
179  mLastError = QObject::tr( "Could not obtain driver for %1" ).arg( mOutputFormat );
180  return CreateOutputError;
181  }
182 
183  gdal::dataset_unique_ptr outputDataset( openOutputFile( outputDriver ) );
184  if ( !outputDataset )
185  {
186  mLastError = QObject::tr( "Could not create output %1" ).arg( mOutputFile );
187  return CreateOutputError;
188  }
189 
190  GDALSetProjection( outputDataset.get(), mOutputCrs.toWkt( QgsCoordinateReferenceSystem::WKT_PREFERRED_GDAL ).toLocal8Bit().data() );
191  GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
192 
193  float outputNodataValue = -FLT_MAX;
194  GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
195 
196 
197  // Take the fast route (process one line at a time) if we can
198  if ( ! requiresMatrix )
199  {
200  // Map of raster names -> blocks
201  std::map<QString, std::unique_ptr<QgsRasterBlock>> inputBlocks;
202  std::map<QString, QgsRasterCalculatorEntry> uniqueRasterEntries;
203  const QList<const QgsRasterCalcNode *> rasterRefNodes = calcNode->findNodes( QgsRasterCalcNode::Type::tRasterRef );
204  for ( const QgsRasterCalcNode *r : rasterRefNodes )
205  {
206  QString layerRef( r->toString().remove( 0, 1 ) );
207  layerRef.chop( 1 );
208  if ( ! inputBlocks.count( layerRef ) )
209  {
210  for ( const QgsRasterCalculatorEntry &ref : std::as_const( mRasterEntries ) )
211  {
212  if ( ref.ref == layerRef )
213  {
214  uniqueRasterEntries[layerRef] = ref;
215  inputBlocks[layerRef ] = std::make_unique<QgsRasterBlock>();
216  }
217  }
218  }
219  }
220 
221  //read / write line by line
222  QMap<QString, QgsRasterBlock * > _rasterData;
223  // Cast to float
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 )
227  {
228  if ( feedback )
229  {
230  feedback->setProgress( 100.0 * static_cast< double >( row ) / mNumOutputRows );
231  }
232 
233  if ( feedback && feedback->isCanceled() )
234  {
235  break;
236  }
237 
238  // Calculates the rect for a single row read
239  QgsRectangle rect( mOutputRectangle );
240  rect.setYMaximum( rect.yMaximum() - rowHeight * row );
241  rect.setYMinimum( rect.yMaximum() - rowHeight );
242 
243  // Read rows into input blocks
244  for ( auto &layerRef : inputBlocks )
245  {
246  QgsRasterCalculatorEntry ref = uniqueRasterEntries[layerRef.first];
247  if ( ref.raster->crs() != mOutputCrs )
248  {
249  QgsRasterProjector proj;
250  proj.setCrs( ref.raster->crs(), mOutputCrs, mTransformContext );
251  proj.setInput( ref.raster->dataProvider() );
253  layerRef.second.reset( proj.block( ref.bandNumber, rect, mNumOutputColumns, 1 ) );
254  }
255  else
256  {
257  layerRef.second.reset( ref.raster->dataProvider()->block( ref.bandNumber, rect, mNumOutputColumns, 1 ) );
258  }
259  }
260 
261  // 1 row X mNumOutputColumns matrix
262  QgsRasterMatrix resultMatrix( mNumOutputColumns, 1, nullptr, outputNodataValue );
263 
264  _rasterData.clear();
265  for ( const auto &layerRef : inputBlocks )
266  {
267  _rasterData.insert( layerRef.first, layerRef.second.get() );
268  }
269 
270  if ( calcNode->calculate( _rasterData, resultMatrix, 0 ) )
271  {
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 )
274  {
275  QgsDebugMsg( QStringLiteral( "RasterIO error!" ) );
276  }
277  }
278  else
279  {
280  //delete the dataset without closing (because it is faster)
281  gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
282  return CalculationError;
283  }
284  }
285 
286  if ( feedback )
287  {
288  feedback->setProgress( 100.0 );
289  }
290  }
291  else // Original code (memory inefficient route)
292  {
293  QMap< QString, QgsRasterBlock * > inputBlocks;
294  QVector<QgsRasterCalculatorEntry>::const_iterator it = mRasterEntries.constBegin();
295  for ( ; it != mRasterEntries.constEnd(); ++it )
296  {
297 
298  std::unique_ptr< QgsRasterBlock > block;
299  // if crs transform needed
300  if ( it->raster->crs() != mOutputCrs )
301  {
302  QgsRasterProjector proj;
303  proj.setCrs( it->raster->crs(), mOutputCrs, it->raster->transformContext() );
304  proj.setInput( it->raster->dataProvider() );
306 
307  QgsRasterBlockFeedback *rasterBlockFeedback = new QgsRasterBlockFeedback();
308  QObject::connect( feedback, &QgsFeedback::canceled, rasterBlockFeedback, &QgsRasterBlockFeedback::cancel );
309  block.reset( proj.block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows, rasterBlockFeedback ) );
310  if ( rasterBlockFeedback->isCanceled() )
311  {
312  qDeleteAll( inputBlocks );
313  return Canceled;
314  }
315  }
316  else
317  {
318  block.reset( it->raster->dataProvider()->block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows ) );
319  }
320  if ( block->isEmpty() )
321  {
322  mLastError = QObject::tr( "Could not allocate required memory for %1" ).arg( it->ref );
323  qDeleteAll( inputBlocks );
324  return MemoryError;
325  }
326  inputBlocks.insert( it->ref, block.release() );
327  }
328 
329  QgsRasterMatrix resultMatrix;
330  resultMatrix.setNodataValue( outputNodataValue );
331 
332  //read / write line by line
333  for ( int i = 0; i < mNumOutputRows; ++i )
334  {
335  if ( feedback )
336  {
337  feedback->setProgress( 100.0 * static_cast< double >( i ) / mNumOutputRows );
338  }
339 
340  if ( feedback && feedback->isCanceled() )
341  {
342  break;
343  }
344 
345  if ( calcNode->calculate( inputBlocks, resultMatrix, i ) )
346  {
347  bool resultIsNumber = resultMatrix.isNumber();
348  float *calcData = new float[mNumOutputColumns];
349 
350  for ( int j = 0; j < mNumOutputColumns; ++j )
351  {
352  calcData[j] = ( float )( resultIsNumber ? resultMatrix.number() : resultMatrix.data()[j] );
353  }
354 
355  //write scanline to the dataset
356  if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
357  {
358  QgsDebugMsg( QStringLiteral( "RasterIO error!" ) );
359  }
360 
361  delete[] calcData;
362  }
363  else
364  {
365  qDeleteAll( inputBlocks );
366  inputBlocks.clear();
367  gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
368  return CalculationError;
369  }
370 
371  }
372 
373  if ( feedback )
374  {
375  feedback->setProgress( 100.0 );
376  }
377 
378  //close datasets and release memory
379  calcNode.reset();
380  qDeleteAll( inputBlocks );
381  inputBlocks.clear();
382 
383  }
384 
385  if ( feedback && feedback->isCanceled() )
386  {
387  //delete the dataset without closing (because it is faster)
388  gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
389  return Canceled;
390  }
391 
392  GDALComputeRasterStatistics( outputRasterBand, true, nullptr, nullptr, nullptr, nullptr, GdalProgressCallback, feedback );
393 
394  return Success;
395 }
396 
397 #ifdef HAVE_OPENCL
398 QgsRasterCalculator::Result QgsRasterCalculator::processCalculationGPU( std::unique_ptr< QgsRasterCalcNode > calcNode, QgsFeedback *feedback )
399 {
400 
401  QString cExpression( calcNode->toString( true ) );
402 
403  QList<const QgsRasterCalcNode *> nodeList( calcNode->findNodes( QgsRasterCalcNode::Type::tRasterRef ) );
404  QSet<QString> capturedTexts;
405  for ( const auto &r : std::as_const( nodeList ) )
406  {
407  QString s( r->toString().remove( 0, 1 ) );
408  s.chop( 1 );
409  capturedTexts.insert( s );
410  }
411 
412  // Extract all references
413  struct LayerRef
414  {
415  QString name;
416  int band;
417  QgsRasterLayer *layer = nullptr;
418  QString varName;
419  QString typeName;
420  size_t index;
421  size_t bufferSize;
422  size_t dataSize;
423  };
424 
425  // Collects all layers, band, name, varName and size information
426  std::vector<LayerRef> inputRefs;
427  size_t refCounter = 0;
428  for ( const auto &r : capturedTexts )
429  {
430  if ( r.startsWith( '"' ) )
431  continue;
432  QStringList parts( r.split( '@' ) );
433  if ( parts.count() != 2 )
434  continue;
435  bool ok = false;
436  LayerRef entry;
437  entry.name = r;
438  entry.band = parts[1].toInt( &ok );
439  for ( const auto &ref : std::as_const( mRasterEntries ) )
440  {
441  if ( ref.ref == entry.name )
442  entry.layer = ref.raster;
443  }
444  if ( !( entry.layer && entry.layer->dataProvider() && ok ) )
445  continue;
446  entry.dataSize = entry.layer->dataProvider()->dataTypeSize( entry.band );
447  switch ( entry.layer->dataProvider()->dataType( entry.band ) )
448  {
450  entry.typeName = QStringLiteral( "unsigned char" );
451  break;
453  entry.typeName = QStringLiteral( "unsigned int" );
454  break;
456  entry.typeName = QStringLiteral( "short" );
457  break;
459  entry.typeName = QStringLiteral( "unsigned int" );
460  break;
462  entry.typeName = QStringLiteral( "int" );
463  break;
465  entry.typeName = QStringLiteral( "float" );
466  break;
467  // FIXME: not sure all OpenCL implementations support double
468  // maybe safer to fall back to the CPU implementation
469  // after a compatibility check
471  entry.typeName = QStringLiteral( "double" );
472  break;
473  default:
474  return BandError;
475  }
476  entry.bufferSize = entry.dataSize * mNumOutputColumns;
477  entry.index = refCounter;
478  entry.varName = QStringLiteral( "input_raster_%1_band_%2" )
479  .arg( refCounter++ )
480  .arg( entry.band );
481  inputRefs.push_back( entry );
482  }
483 
484  // May throw an openCL exception
485  try
486  {
487  // Prepare context and queue
488  cl::Context ctx( QgsOpenClUtils::context() );
489  cl::CommandQueue queue( QgsOpenClUtils::commandQueue() );
490 
491  // Create the C expression
492  std::vector<cl::Buffer> inputBuffers;
493  inputBuffers.reserve( inputRefs.size() );
494  QStringList inputArgs;
495  for ( const auto &ref : inputRefs )
496  {
497  cExpression.replace( QStringLiteral( "\"%1\"" ).arg( ref.name ), QStringLiteral( "%1[i]" ).arg( ref.varName ) );
498  inputArgs.append( QStringLiteral( "__global %1 *%2" )
499  .arg( ref.typeName, ref.varName ) );
500  inputBuffers.push_back( cl::Buffer( ctx, CL_MEM_READ_ONLY, ref.bufferSize, nullptr, nullptr ) );
501  }
502 
503  //qDebug() << cExpression;
504 
505  // Create the program
506  QString programTemplate( R"CL(
507  // Inputs:
508  ##INPUT_DESC##
509  // Expression: ##EXPRESSION_ORIGINAL##
510  __kernel void rasterCalculator( ##INPUT##
511  __global float *resultLine
512  )
513  {
514  // Get the index of the current element
515  const int i = get_global_id(0);
516  // Check for nodata in input
517  if ( ##INPUT_NODATA_CHECK## )
518  resultLine[i] = -FLT_MAX;
519  // Expression
520  else
521  resultLine[i] = ##EXPRESSION##;
522  }
523  )CL" );
524 
525  QStringList inputDesc;
526  QStringList inputNoDataCheck;
527  for ( const auto &ref : inputRefs )
528  {
529  inputDesc.append( QStringLiteral( " // %1 = %2" ).arg( ref.varName, ref.name ) );
530  if ( ref.layer->dataProvider()->sourceHasNoDataValue( ref.band ) )
531  {
532  inputNoDataCheck.append( QStringLiteral( "(float) %1[i] == (float) %2" ).arg( ref.varName, QString::number( ref.layer->dataProvider()->sourceNoDataValue( ref.band ), 'g', 10 ) ) );
533  }
534  }
535 
536  programTemplate = programTemplate.replace( QLatin1String( "##INPUT_NODATA_CHECK##" ), inputNoDataCheck.isEmpty() ? QStringLiteral( "false" ) : inputNoDataCheck.join( QLatin1String( " || " ) ) );
537  programTemplate = programTemplate.replace( QLatin1String( "##INPUT_DESC##" ), inputDesc.join( '\n' ) );
538  programTemplate = programTemplate.replace( QLatin1String( "##INPUT##" ), !inputArgs.isEmpty() ? ( inputArgs.join( ',' ).append( ',' ) ) : QChar( ' ' ) );
539  programTemplate = programTemplate.replace( QLatin1String( "##EXPRESSION##" ), cExpression );
540  programTemplate = programTemplate.replace( QLatin1String( "##EXPRESSION_ORIGINAL##" ), calcNode->toString( ) );
541 
542  //qDebug() << programTemplate;
543 
544  // Create a program from the kernel source
545  cl::Program program( QgsOpenClUtils::buildProgram( programTemplate, QgsOpenClUtils::ExceptionBehavior::Throw ) );
546 
547  // Create the buffers, output is float32 (4 bytes)
548  // We assume size of float = 4 because that's the size used by OpenCL and IEEE 754
549  Q_ASSERT( sizeof( float ) == 4 );
550  std::size_t resultBufferSize( 4 * static_cast<size_t>( mNumOutputColumns ) );
551  cl::Buffer resultLineBuffer( ctx, CL_MEM_WRITE_ONLY,
552  resultBufferSize, nullptr, nullptr );
553 
554  auto kernel = cl::Kernel( program, "rasterCalculator" );
555 
556  for ( unsigned int i = 0; i < inputBuffers.size() ; i++ )
557  {
558  kernel.setArg( i, inputBuffers.at( i ) );
559  }
560  kernel.setArg( static_cast<unsigned int>( inputBuffers.size() ), resultLineBuffer );
561 
562  QgsOpenClUtils::CPLAllocator<float> resultLine( static_cast<size_t>( mNumOutputColumns ) );
563 
564  //open output dataset for writing
565  GDALDriverH outputDriver = openOutputDriver();
566  if ( !outputDriver )
567  {
568  mLastError = QObject::tr( "Could not obtain driver for %1" ).arg( mOutputFormat );
569  return CreateOutputError;
570  }
571 
572  gdal::dataset_unique_ptr outputDataset( openOutputFile( outputDriver ) );
573  if ( !outputDataset )
574  {
575  mLastError = QObject::tr( "Could not create output %1" ).arg( mOutputFile );
576  return CreateOutputError;
577  }
578 
579  GDALSetProjection( outputDataset.get(), mOutputCrs.toWkt( QgsCoordinateReferenceSystem::WKT_PREFERRED_GDAL ).toLocal8Bit().data() );
580 
581 
582  GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
583  if ( !outputRasterBand )
584  return BandError;
585 
586  const float outputNodataValue = -FLT_MAX;
587  GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
588 
589  // Input block (buffer)
590  std::unique_ptr<QgsRasterBlock> block;
591 
592  // Run kernel on all scanlines
593  auto rowHeight = mOutputRectangle.height() / mNumOutputRows;
594  for ( int line = 0; line < mNumOutputRows; line++ )
595  {
596  if ( feedback && feedback->isCanceled() )
597  {
598  break;
599  }
600 
601  if ( feedback )
602  {
603  feedback->setProgress( 100.0 * static_cast< double >( line ) / mNumOutputRows );
604  }
605 
606  // Read lines from rasters into the buffers
607  for ( const auto &ref : inputRefs )
608  {
609  // Read one row
610  QgsRectangle rect( mOutputRectangle );
611  rect.setYMaximum( rect.yMaximum() - rowHeight * line );
612  rect.setYMinimum( rect.yMaximum() - rowHeight );
613 
614  // TODO: check if this is too slow
615  // if crs transform needed
616  if ( ref.layer->crs() != mOutputCrs )
617  {
618  QgsRasterProjector proj;
619  proj.setCrs( ref.layer->crs(), mOutputCrs, ref.layer->transformContext() );
620  proj.setInput( ref.layer->dataProvider() );
622  block.reset( proj.block( ref.band, rect, mNumOutputColumns, 1 ) );
623  }
624  else
625  {
626  block.reset( ref.layer->dataProvider()->block( ref.band, rect, mNumOutputColumns, 1 ) );
627  }
628 
629  //for ( int i = 0; i < mNumOutputColumns; i++ )
630  // qDebug() << "Input: " << line << i << ref.varName << " = " << block->value( 0, i );
631  //qDebug() << "Writing buffer " << ref.index;
632 
633  Q_ASSERT( ref.bufferSize == static_cast<size_t>( block->data().size( ) ) );
634  queue.enqueueWriteBuffer( inputBuffers[ref.index], CL_TRUE, 0,
635  ref.bufferSize, block->bits() );
636 
637  }
638  // Run the kernel
639  queue.enqueueNDRangeKernel(
640  kernel,
641  0,
642  cl::NDRange( mNumOutputColumns )
643  );
644 
645  // Write the result
646  queue.enqueueReadBuffer( resultLineBuffer, CL_TRUE, 0,
647  resultBufferSize, resultLine.get() );
648 
649  //for ( int i = 0; i < mNumOutputColumns; i++ )
650  // qDebug() << "Output: " << line << i << " = " << resultLine[i];
651 
652  if ( GDALRasterIO( outputRasterBand, GF_Write, 0, line, mNumOutputColumns, 1, resultLine.get(), mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
653  {
654  return CreateOutputError;
655  }
656  }
657 
658  if ( feedback && feedback->isCanceled() )
659  {
660  //delete the dataset without closing (because it is faster)
661  gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
662  return Canceled;
663  }
664 
665  inputBuffers.clear();
666 
667  GDALComputeRasterStatistics( outputRasterBand, true, nullptr, nullptr, nullptr, nullptr, GdalProgressCallback, feedback );
668 
669  }
670  catch ( cl::Error &e )
671  {
672  mLastError = e.what();
673  return CreateOutputError;
674  }
675 
676  return Success;
677 }
678 #endif
679 
680 GDALDriverH QgsRasterCalculator::openOutputDriver()
681 {
682  //open driver
683  GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
684 
685  if ( !outputDriver )
686  {
687  return outputDriver; //return nullptr, driver does not exist
688  }
689 
690  if ( !QgsGdalUtils::supportsRasterCreate( outputDriver ) )
691  {
692  return nullptr; //driver exist, but it does not support the create operation
693  }
694 
695  return outputDriver;
696 }
697 
698 gdal::dataset_unique_ptr QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver )
699 {
700  //open output file
701  char **papszOptions = nullptr;
702  gdal::dataset_unique_ptr outputDataset( GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions ) );
703  if ( !outputDataset )
704  {
705  return nullptr;
706  }
707 
708  //assign georef information
709  double geotransform[6];
710  outputGeoTransform( geotransform );
711  GDALSetGeoTransform( outputDataset.get(), geotransform );
712 
713  return outputDataset;
714 }
715 
716 void QgsRasterCalculator::outputGeoTransform( double *transform ) const
717 {
718  transform[0] = mOutputRectangle.xMinimum();
719  transform[1] = mOutputRectangle.width() / mNumOutputColumns;
720  transform[2] = 0;
721  transform[3] = mOutputRectangle.yMaximum();
722  transform[4] = 0;
723  transform[5] = -mOutputRectangle.height() / mNumOutputRows;
724 }
725 
727 {
728  return mLastError;
729 }
730 
731 QVector<QgsRasterCalculatorEntry> QgsRasterCalculatorEntry::rasterEntries()
732 {
733  QVector<QgsRasterCalculatorEntry> availableEntries;
734  const QMap<QString, QgsMapLayer *> &layers = QgsProject::instance()->mapLayers();
735 
736  auto uniqueRasterBandIdentifier = [ & ]( QgsRasterCalculatorEntry & entry ) -> bool
737  {
738  unsigned int i( 1 );
739  entry.ref = QStringLiteral( "%[email protected]%2" ).arg( entry.raster->name() ).arg( entry.bandNumber );
740  while ( true )
741  {
742  bool unique( true );
743  for ( const auto &ref : std::as_const( availableEntries ) )
744  {
745  // Safety belt
746  if ( !( entry.raster && ref.raster ) )
747  continue;
748  // Check if is another band of the same raster
749  if ( ref.raster->publicSource() == entry.raster->publicSource() )
750  {
751  if ( ref.bandNumber != entry.bandNumber )
752  {
753  continue;
754  }
755  else // a layer with the same data source was already added to the list
756  {
757  return false;
758  }
759  }
760  // If same name but different source
761  if ( ref.ref == entry.ref )
762  {
763  unique = false;
764  entry.ref = QStringLiteral( "%1_%[email protected]%3" ).arg( entry.raster->name() ).arg( i++ ).arg( entry.bandNumber );
765  }
766  }
767  if ( unique )
768  return true;
769  }
770  };
771 
772  QMap<QString, QgsMapLayer *>::const_iterator layerIt = layers.constBegin();
773  for ( ; layerIt != layers.constEnd(); ++layerIt )
774  {
775  QgsRasterLayer *rlayer = qobject_cast<QgsRasterLayer *>( layerIt.value() );
776  if ( rlayer && rlayer->dataProvider() && ( rlayer->dataProvider()->capabilities() & QgsRasterDataProvider::Size ) )
777  {
778  //get number of bands
779  for ( int i = 0; i < rlayer->bandCount(); ++i )
780  {
782  entry.raster = rlayer;
783  entry.bandNumber = i + 1;
784  if ( ! uniqueRasterBandIdentifier( entry ) )
785  break;
786  availableEntries.push_back( entry );
787  }
788  }
789  }
790  return availableEntries;
791 }
@ Float32
Thirty two bit floating point (float)
@ Int16
Sixteen bit signed integer (qint16)
@ UInt16
Sixteen bit unsigned integer (quint16)
@ Byte
Eight bit unsigned integer (quint8)
@ Int32
Thirty two bit signed integer (qint32)
@ Float64
Sixty four bit floating point (double)
@ UInt32
Thirty two bit unsigned integer (quint32)
This class represents a coordinate reference system (CRS).
@ WKT_PREFERRED_GDAL
Preferred format for conversion of CRS to WKT for use with the GDAL library.
QString toWkt(WktVariant variant=WKT1_GDAL, 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.
Definition: qgsfeedback.h:45
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 ...
Definition: qgsfeedback.h:108
bool isCanceled() const SIP_HOLDGIL
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
static bool supportsRasterCreate(GDALDriverH driver)
Reads whether a driver supports GDALCreate() for raster purposes.
QgsCoordinateReferenceSystem crs
Definition: qgsmaplayer.h:79
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.
Encapsulates a QGIS project, including sets of map layers and their styles, layouts,...
Definition: qgsproject.h:104
static QgsProject * instance()
Returns the QgsProject singleton instance.
Definition: qgsproject.cpp:479
QgsCoordinateTransformContext transformContext
Definition: qgsproject.h:110
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.
@ Size
Original data source size (and thus resolution) is known, it is not always available,...
virtual int capabilities() const
Returns a bitmask containing the supported capabilities.
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.
double * data()
Returns data array (but not ownership)
bool isNumber() const
Returns true if matrix is 1x1 (=scalar number)
void setNodataValue(double d)
double number() const
QgsRasterProjector implements approximate projection support for it calculates grid of points in sour...
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.
Definition: qgsrectangle.h:42
double yMaximum() const SIP_HOLDGIL
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:193
double xMinimum() const SIP_HOLDGIL
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:188
void setYMinimum(double y) SIP_HOLDGIL
Set the minimum y value.
Definition: qgsrectangle.h:161
double height() const SIP_HOLDGIL
Returns the height of the rectangle.
Definition: qgsrectangle.h:230
void setYMaximum(double y) SIP_HOLDGIL
Set the maximum y value.
Definition: qgsrectangle.h:166
double width() const SIP_HOLDGIL
Returns the width of the rectangle.
Definition: qgsrectangle.h:223
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.
Definition: qgsogrutils.cpp:80
std::unique_ptr< std::remove_pointer< GDALDatasetH >::type, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
Definition: qgsogrutils.h:139
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
int CPL_STDCALL GdalProgressCallback(double dfComplete, const char *pszMessage, void *pProgressArg)
const QgsCoordinateReferenceSystem & outputCrs
const QString & typeName
Tiny smart-pointer-like wrapper around CPLMalloc and CPLFree: this is needed because OpenCL C++ API m...