QGIS API Documentation  3.14.0-Pi (9f7028fd23)
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 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 )
40  : mFormulaString( formulaString )
41  , mOutputFile( outputFile )
42  , mOutputFormat( outputFormat )
43  , mOutputRectangle( outputExtent )
44  , mNumOutputColumns( nOutputColumns )
45  , mNumOutputRows( nOutputRows )
46  , mRasterEntries( rasterEntries )
47  , mTransformContext( transformContext )
48 {
49 
50 }
51 
52 QgsRasterCalculator::QgsRasterCalculator( const QString &formulaString, const QString &outputFile, const QString &outputFormat,
53  const QgsRectangle &outputExtent, const QgsCoordinateReferenceSystem &outputCrs, int nOutputColumns, int nOutputRows,
54  const QVector<QgsRasterCalculatorEntry> &rasterEntries, const QgsCoordinateTransformContext &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 )
64 {
65 
66 }
67 
68 // Deprecated!
69 QgsRasterCalculator::QgsRasterCalculator( const QString &formulaString, const QString &outputFile, const QString &outputFormat,
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 )
78 {
79  //default to first layer's crs
80  mOutputCrs = mRasterEntries.at( 0 ).raster->crs();
81  mTransformContext = QgsProject::instance()->transformContext();
82 }
83 
84 
85 // Deprecated!
86 QgsRasterCalculator::QgsRasterCalculator( const QString &formulaString, const QString &outputFile, const QString &outputFormat,
87  const QgsRectangle &outputExtent, const QgsCoordinateReferenceSystem &outputCrs, int nOutputColumns, int nOutputRows, const QVector<QgsRasterCalculatorEntry> &rasterEntries )
88  : mFormulaString( formulaString )
89  , mOutputFile( outputFile )
90  , mOutputFormat( outputFormat )
91  , mOutputRectangle( outputExtent )
92  , mOutputCrs( outputCrs )
93  , mNumOutputColumns( nOutputColumns )
94  , mNumOutputRows( nOutputRows )
95  , mRasterEntries( rasterEntries )
96 {
97  mTransformContext = QgsProject::instance()->transformContext();
98 }
99 
101 {
102  mLastError.clear();
103 
104  //prepare search string / tree
105  std::unique_ptr< QgsRasterCalcNode > calcNode( QgsRasterCalcNode::parseRasterCalcString( mFormulaString, mLastError ) );
106  if ( !calcNode )
107  {
108  //error
109  return ParserError;
110  }
111 
112  // Check input layers and bands
113  for ( const auto &entry : qgis::as_const( mRasterEntries ) )
114  {
115  if ( !entry.raster ) // no raster layer in entry
116  {
117  mLastError = QObject::tr( "No raster layer for entry %1" ).arg( entry.ref );
118  return InputLayerError;
119  }
120  if ( entry.bandNumber <= 0 || entry.bandNumber > entry.raster->bandCount() )
121  {
122  mLastError = QObject::tr( "Band number %1 is not valid for entry %2" ).arg( entry.bandNumber ).arg( entry.ref );
123  return BandError;
124  }
125  }
126 
127  // Check if we need to read the raster as a whole (which is memory inefficient
128  // and not interruptible by the user) by checking if any raster matrix nodes are
129  // in the expression
130  bool requiresMatrix = ! calcNode->findNodes( QgsRasterCalcNode::Type::tMatrix ).isEmpty();
131 
132 #ifdef HAVE_OPENCL
133  // Check for matrix nodes, GPU implementation does not support them
134  QList<const QgsRasterCalcNode *> nodeList;
135  if ( QgsOpenClUtils::enabled() && QgsOpenClUtils::available() && ! requiresMatrix )
136  {
137  return processCalculationGPU( std::move( calcNode ), feedback );
138  }
139 #endif
140 
141  //open output dataset for writing
142  GDALDriverH outputDriver = openOutputDriver();
143  if ( !outputDriver )
144  {
145  mLastError = QObject::tr( "Could not obtain driver for %1" ).arg( mOutputFormat );
146  return CreateOutputError;
147  }
148 
149  gdal::dataset_unique_ptr outputDataset( openOutputFile( outputDriver ) );
150  if ( !outputDataset )
151  {
152  mLastError = QObject::tr( "Could not create output %1" ).arg( mOutputFile );
153  return CreateOutputError;
154  }
155 
156  GDALSetProjection( outputDataset.get(), mOutputCrs.toWkt( QgsCoordinateReferenceSystem::WKT_PREFERRED_GDAL ).toLocal8Bit().data() );
157  GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
158 
159  float outputNodataValue = -FLT_MAX;
160  GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
161 
162 
163  // Take the fast route (process one line at a time) if we can
164  if ( ! requiresMatrix )
165  {
166  // Map of raster names -> blocks
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 );
170  for ( const QgsRasterCalcNode *r : rasterRefNodes )
171  {
172  QString layerRef( r->toString().remove( 0, 1 ) );
173  layerRef.chop( 1 );
174  if ( ! inputBlocks.count( layerRef ) )
175  {
176  for ( const QgsRasterCalculatorEntry &ref : qgis::as_const( mRasterEntries ) )
177  {
178  if ( ref.ref == layerRef )
179  {
180  uniqueRasterEntries[layerRef] = ref;
181  inputBlocks[layerRef ] = qgis::make_unique<QgsRasterBlock>();
182  }
183  }
184  }
185  }
186 
187  //read / write line by line
188  QMap<QString, QgsRasterBlock * > _rasterData;
189  // Cast to float
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 )
193  {
194  if ( feedback )
195  {
196  feedback->setProgress( 100.0 * static_cast< double >( row ) / mNumOutputRows );
197  }
198 
199  if ( feedback && feedback->isCanceled() )
200  {
201  break;
202  }
203 
204  // Calculates the rect for a single row read
205  QgsRectangle rect( mOutputRectangle );
206  rect.setYMaximum( rect.yMaximum() - rowHeight * row );
207  rect.setYMinimum( rect.yMaximum() - rowHeight );
208 
209  // Read rows into input blocks
210  for ( auto &layerRef : inputBlocks )
211  {
212  QgsRasterCalculatorEntry ref = uniqueRasterEntries[layerRef.first];
213  if ( ref.raster->crs() != mOutputCrs )
214  {
215  QgsRasterProjector proj;
216  proj.setCrs( ref.raster->crs(), mOutputCrs, mTransformContext );
217  proj.setInput( ref.raster->dataProvider() );
219  layerRef.second.reset( proj.block( ref.bandNumber, rect, mNumOutputColumns, 1 ) );
220  }
221  else
222  {
223  layerRef.second.reset( ref.raster->dataProvider()->block( ref.bandNumber, rect, mNumOutputColumns, 1 ) );
224  }
225  }
226 
227  // 1 row X mNumOutputColumns matrix
228  QgsRasterMatrix resultMatrix( mNumOutputColumns, 1, nullptr, outputNodataValue );
229 
230  _rasterData.clear();
231  for ( const auto &layerRef : inputBlocks )
232  {
233  _rasterData.insert( layerRef.first, layerRef.second.get() );
234  }
235 
236  if ( calcNode->calculate( _rasterData, resultMatrix, 0 ) )
237  {
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 )
240  {
241  QgsDebugMsg( QStringLiteral( "RasterIO error!" ) );
242  }
243  }
244  else
245  {
246  //delete the dataset without closing (because it is faster)
247  gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
248  return CalculationError;
249  }
250  }
251 
252  if ( feedback )
253  {
254  feedback->setProgress( 100.0 );
255  }
256  }
257  else // Original code (memory inefficient route)
258  {
259  QMap< QString, QgsRasterBlock * > inputBlocks;
260  QVector<QgsRasterCalculatorEntry>::const_iterator it = mRasterEntries.constBegin();
261  for ( ; it != mRasterEntries.constEnd(); ++it )
262  {
263 
264  std::unique_ptr< QgsRasterBlock > block;
265  // if crs transform needed
266  if ( it->raster->crs() != mOutputCrs )
267  {
268  QgsRasterProjector proj;
269  proj.setCrs( it->raster->crs(), mOutputCrs, it->raster->transformContext() );
270  proj.setInput( it->raster->dataProvider() );
272 
273  QgsRasterBlockFeedback *rasterBlockFeedback = new QgsRasterBlockFeedback();
274  QObject::connect( feedback, &QgsFeedback::canceled, rasterBlockFeedback, &QgsRasterBlockFeedback::cancel );
275  block.reset( proj.block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows, rasterBlockFeedback ) );
276  if ( rasterBlockFeedback->isCanceled() )
277  {
278  qDeleteAll( inputBlocks );
279  return Canceled;
280  }
281  }
282  else
283  {
284  block.reset( it->raster->dataProvider()->block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows ) );
285  }
286  if ( block->isEmpty() )
287  {
288  mLastError = QObject::tr( "Could not allocate required memory for %1" ).arg( it->ref );
289  qDeleteAll( inputBlocks );
290  return MemoryError;
291  }
292  inputBlocks.insert( it->ref, block.release() );
293  }
294 
295  QgsRasterMatrix resultMatrix;
296  resultMatrix.setNodataValue( outputNodataValue );
297 
298  //read / write line by line
299  for ( int i = 0; i < mNumOutputRows; ++i )
300  {
301  if ( feedback )
302  {
303  feedback->setProgress( 100.0 * static_cast< double >( i ) / mNumOutputRows );
304  }
305 
306  if ( feedback && feedback->isCanceled() )
307  {
308  break;
309  }
310 
311  if ( calcNode->calculate( inputBlocks, resultMatrix, i ) )
312  {
313  bool resultIsNumber = resultMatrix.isNumber();
314  float *calcData = new float[mNumOutputColumns];
315 
316  for ( int j = 0; j < mNumOutputColumns; ++j )
317  {
318  calcData[j] = ( float )( resultIsNumber ? resultMatrix.number() : resultMatrix.data()[j] );
319  }
320 
321  //write scanline to the dataset
322  if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
323  {
324  QgsDebugMsg( QStringLiteral( "RasterIO error!" ) );
325  }
326 
327  delete[] calcData;
328  }
329  else
330  {
331  qDeleteAll( inputBlocks );
332  inputBlocks.clear();
333  gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
334  return CalculationError;
335  }
336 
337  }
338 
339  if ( feedback )
340  {
341  feedback->setProgress( 100.0 );
342  }
343 
344  //close datasets and release memory
345  calcNode.reset();
346  qDeleteAll( inputBlocks );
347  inputBlocks.clear();
348 
349  }
350 
351  if ( feedback && feedback->isCanceled() )
352  {
353  //delete the dataset without closing (because it is faster)
354  gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
355  return Canceled;
356  }
357  return Success;
358 }
359 
360 #ifdef HAVE_OPENCL
361 QgsRasterCalculator::Result QgsRasterCalculator::processCalculationGPU( std::unique_ptr< QgsRasterCalcNode > calcNode, QgsFeedback *feedback )
362 {
363 
364  QString cExpression( calcNode->toString( true ) );
365 
366  QList<const QgsRasterCalcNode *> nodeList( calcNode->findNodes( QgsRasterCalcNode::Type::tRasterRef ) );
367  QSet<QString> capturedTexts;
368  for ( const auto &r : qgis::as_const( nodeList ) )
369  {
370  QString s( r->toString().remove( 0, 1 ) );
371  s.chop( 1 );
372  capturedTexts.insert( s );
373  }
374 
375  // Extract all references
376  struct LayerRef
377  {
378  QString name;
379  int band;
380  QgsRasterLayer *layer = nullptr;
381  QString varName;
382  QString typeName;
383  size_t index;
384  size_t bufferSize;
385  size_t dataSize;
386  };
387 
388  // Collects all layers, band, name, varName and size information
389  std::vector<LayerRef> inputRefs;
390  size_t refCounter = 0;
391  for ( const auto &r : capturedTexts )
392  {
393  if ( r.startsWith( '"' ) )
394  continue;
395  QStringList parts( r.split( '@' ) );
396  if ( parts.count() != 2 )
397  continue;
398  bool ok = false;
399  LayerRef entry;
400  entry.name = r;
401  entry.band = parts[1].toInt( &ok );
402  for ( const auto &ref : mRasterEntries )
403  {
404  if ( ref.ref == entry.name )
405  entry.layer = ref.raster;
406  }
407  if ( !( entry.layer && entry.layer->dataProvider() && ok ) )
408  continue;
409  entry.dataSize = entry.layer->dataProvider()->dataTypeSize( entry.band );
410  switch ( entry.layer->dataProvider()->dataType( entry.band ) )
411  {
412  case Qgis::DataType::Byte:
413  entry.typeName = QStringLiteral( "unsigned char" );
414  break;
415  case Qgis::DataType::UInt16:
416  entry.typeName = QStringLiteral( "unsigned int" );
417  break;
418  case Qgis::DataType::Int16:
419  entry.typeName = QStringLiteral( "short" );
420  break;
421  case Qgis::DataType::UInt32:
422  entry.typeName = QStringLiteral( "unsigned int" );
423  break;
424  case Qgis::DataType::Int32:
425  entry.typeName = QStringLiteral( "int" );
426  break;
427  case Qgis::DataType::Float32:
428  entry.typeName = QStringLiteral( "float" );
429  break;
430  // FIXME: not sure all OpenCL implementations support double
431  // maybe safer to fall back to the CPU implementation
432  // after a compatibility check
433  case Qgis::DataType::Float64:
434  entry.typeName = QStringLiteral( "double" );
435  break;
436  default:
437  return BandError;
438  }
439  entry.bufferSize = entry.dataSize * mNumOutputColumns;
440  entry.index = refCounter;
441  entry.varName = QStringLiteral( "input_raster_%1_band_%2" )
442  .arg( refCounter++ )
443  .arg( entry.band );
444  inputRefs.push_back( entry );
445  }
446 
447  // May throw an openCL exception
448  try
449  {
450  // Prepare context and queue
451  cl::Context ctx( QgsOpenClUtils::context() );
452  cl::CommandQueue queue( QgsOpenClUtils::commandQueue() );
453 
454  // Create the C expression
455  std::vector<cl::Buffer> inputBuffers;
456  inputBuffers.reserve( inputRefs.size() );
457  QStringList inputArgs;
458  for ( const auto &ref : inputRefs )
459  {
460  cExpression.replace( QStringLiteral( "\"%1\"" ).arg( ref.name ), QStringLiteral( "%1[i]" ).arg( ref.varName ) );
461  inputArgs.append( QStringLiteral( "__global %1 *%2" )
462  .arg( ref.typeName )
463  .arg( ref.varName ) );
464  inputBuffers.push_back( cl::Buffer( ctx, CL_MEM_READ_ONLY, ref.bufferSize, nullptr, nullptr ) );
465  }
466 
467  //qDebug() << cExpression;
468 
469  // Create the program
470  QString programTemplate( R"CL(
471  // Inputs:
472  ##INPUT_DESC##
473  // Expression: ##EXPRESSION_ORIGINAL##
474  __kernel void rasterCalculator( ##INPUT##
475  __global float *resultLine
476  )
477  {
478  // Get the index of the current element
479  const int i = get_global_id(0);
480  // Expression
481  resultLine[i] = ##EXPRESSION##;
482  }
483  )CL" );
484 
485  QStringList inputDesc;
486  for ( const auto &ref : inputRefs )
487  {
488  inputDesc.append( QStringLiteral( " // %1 = %2" ).arg( ref.varName ).arg( ref.name ) );
489  }
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( ) );
494 
495  // qDebug() << programTemplate;
496 
497  // Create a program from the kernel source
498  cl::Program program( QgsOpenClUtils::buildProgram( programTemplate, QgsOpenClUtils::ExceptionBehavior::Throw ) );
499 
500  // Create the buffers, output is float32 (4 bytes)
501  // We assume size of float = 4 because that's the size used by OpenCL and IEEE 754
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 );
506 
507  auto kernel = cl::Kernel( program, "rasterCalculator" );
508 
509  for ( unsigned int i = 0; i < inputBuffers.size() ; i++ )
510  {
511  kernel.setArg( i, inputBuffers.at( i ) );
512  }
513  kernel.setArg( static_cast<unsigned int>( inputBuffers.size() ), resultLineBuffer );
514 
515  QgsOpenClUtils::CPLAllocator<float> resultLine( static_cast<size_t>( mNumOutputColumns ) );
516 
517  //open output dataset for writing
518  GDALDriverH outputDriver = openOutputDriver();
519  if ( !outputDriver )
520  {
521  mLastError = QObject::tr( "Could not obtain driver for %1" ).arg( mOutputFormat );
522  return CreateOutputError;
523  }
524 
525  gdal::dataset_unique_ptr outputDataset( openOutputFile( outputDriver ) );
526  if ( !outputDataset )
527  {
528  mLastError = QObject::tr( "Could not create output %1" ).arg( mOutputFile );
529  return CreateOutputError;
530  }
531 
532  GDALSetProjection( outputDataset.get(), mOutputCrs.toWkt( QgsCoordinateReferenceSystem::WKT_PREFERRED_GDAL ).toLocal8Bit().data() );
533 
534 
535  GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
536  if ( !outputRasterBand )
537  return BandError;
538 
539  float outputNodataValue = -FLT_MAX;
540  GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
541 
542  // Input block (buffer)
543  std::unique_ptr<QgsRasterBlock> block;
544 
545  // Run kernel on all scanlines
546  auto rowHeight = mOutputRectangle.height() / mNumOutputRows;
547  for ( int line = 0; line < mNumOutputRows; line++ )
548  {
549  if ( feedback && feedback->isCanceled() )
550  {
551  break;
552  }
553 
554  if ( feedback )
555  {
556  feedback->setProgress( 100.0 * static_cast< double >( line ) / mNumOutputRows );
557  }
558 
559  // Read lines from rasters into the buffers
560  for ( const auto &ref : inputRefs )
561  {
562  // Read one row
563  QgsRectangle rect( mOutputRectangle );
564  rect.setYMaximum( rect.yMaximum() - rowHeight * line );
565  rect.setYMinimum( rect.yMaximum() - rowHeight );
566 
567  // TODO: check if this is too slow
568  // if crs transform needed
569  if ( ref.layer->crs() != mOutputCrs )
570  {
571  QgsRasterProjector proj;
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 ) );
576  }
577  else
578  {
579  block.reset( ref.layer->dataProvider()->block( ref.band, rect, mNumOutputColumns, 1 ) );
580  }
581 
582  //for ( int i = 0; i < mNumOutputColumns; i++ )
583  // qDebug() << "Input: " << line << i << ref.varName << " = " << block->value( 0, i );
584  //qDebug() << "Writing buffer " << ref.index;
585 
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() );
589 
590  }
591  // Run the kernel
592  queue.enqueueNDRangeKernel(
593  kernel,
594  0,
595  cl::NDRange( mNumOutputColumns )
596  );
597 
598  // Write the result
599  queue.enqueueReadBuffer( resultLineBuffer, CL_TRUE, 0,
600  resultBufferSize, resultLine.get() );
601 
602  //for ( int i = 0; i < mNumOutputColumns; i++ )
603  // qDebug() << "Output: " << line << i << " = " << resultLine[i];
604 
605  if ( GDALRasterIO( outputRasterBand, GF_Write, 0, line, mNumOutputColumns, 1, resultLine.get(), mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
606  {
607  return CreateOutputError;
608  }
609  }
610 
611  if ( feedback && feedback->isCanceled() )
612  {
613  //delete the dataset without closing (because it is faster)
614  gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
615  return Canceled;
616  }
617 
618  inputBuffers.clear();
619 
620  }
621  catch ( cl::Error &e )
622  {
623  mLastError = e.what();
624  return CreateOutputError;
625  }
626 
627  return Success;
628 }
629 #endif
630 
631 GDALDriverH QgsRasterCalculator::openOutputDriver()
632 {
633  //open driver
634  GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
635 
636  if ( !outputDriver )
637  {
638  return outputDriver; //return nullptr, driver does not exist
639  }
640 
641  if ( !QgsGdalUtils::supportsRasterCreate( outputDriver ) )
642  {
643  return nullptr; //driver exist, but it does not support the create operation
644  }
645 
646  return outputDriver;
647 }
648 
649 gdal::dataset_unique_ptr QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver )
650 {
651  //open output file
652  char **papszOptions = nullptr;
653  gdal::dataset_unique_ptr outputDataset( GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions ) );
654  if ( !outputDataset )
655  {
656  return nullptr;
657  }
658 
659  //assign georef information
660  double geotransform[6];
661  outputGeoTransform( geotransform );
662  GDALSetGeoTransform( outputDataset.get(), geotransform );
663 
664  return outputDataset;
665 }
666 
667 void QgsRasterCalculator::outputGeoTransform( double *transform ) const
668 {
669  transform[0] = mOutputRectangle.xMinimum();
670  transform[1] = mOutputRectangle.width() / mNumOutputColumns;
671  transform[2] = 0;
672  transform[3] = mOutputRectangle.yMaximum();
673  transform[4] = 0;
674  transform[5] = -mOutputRectangle.height() / mNumOutputRows;
675 }
676 
678 {
679  return mLastError;
680 }
681 
682 QVector<QgsRasterCalculatorEntry> QgsRasterCalculatorEntry::rasterEntries()
683 {
684  QVector<QgsRasterCalculatorEntry> availableEntries;
685  const QMap<QString, QgsMapLayer *> &layers = QgsProject::instance()->mapLayers();
686 
687  auto uniqueRasterBandIdentifier = [ & ]( QgsRasterCalculatorEntry & entry ) -> bool
688  {
689  unsigned int i( 1 );
690  entry.ref = QStringLiteral( "%1@%2" ).arg( entry.raster->name() ).arg( entry.bandNumber );
691  while ( true )
692  {
693  bool unique( true );
694  for ( const auto &ref : qgis::as_const( availableEntries ) )
695  {
696  // Safety belt
697  if ( !( entry.raster && ref.raster ) )
698  continue;
699  // Check if is another band of the same raster
700  if ( ref.raster->publicSource() == entry.raster->publicSource() )
701  {
702  if ( ref.bandNumber != entry.bandNumber )
703  {
704  continue;
705  }
706  else // a layer with the same data source was already added to the list
707  {
708  return false;
709  }
710  }
711  // If same name but different source
712  if ( ref.ref == entry.ref )
713  {
714  unique = false;
715  entry.ref = QStringLiteral( "%1_%2@%3" ).arg( entry.raster->name() ).arg( i++ ).arg( entry.bandNumber );
716  }
717  }
718  if ( unique )
719  return true;
720  }
721  };
722 
723  QMap<QString, QgsMapLayer *>::const_iterator layerIt = layers.constBegin();
724  for ( ; layerIt != layers.constEnd(); ++layerIt )
725  {
726  QgsRasterLayer *rlayer = qobject_cast<QgsRasterLayer *>( layerIt.value() );
727  if ( rlayer && rlayer->dataProvider() && rlayer->providerType() == QLatin1String( "gdal" ) )
728  {
729  //get number of bands
730  for ( int i = 0; i < rlayer->bandCount(); ++i )
731  {
733  entry.raster = rlayer;
734  entry.bandNumber = i + 1;
735  if ( ! uniqueRasterBandIdentifier( entry ) )
736  break;
737  availableEntries.push_back( entry );
738  }
739  }
740  }
741  return availableEntries;
742 }
QgsMapLayer::crs
QgsCoordinateReferenceSystem crs
Definition: qgsmaplayer.h:88
QgsRasterCalculatorEntry::raster
QgsRasterLayer * raster
Raster layer associated with entry.
Definition: qgsrastercalculator.h:64
QgsFeedback::setProgress
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:75
qgsrasterprojector.h
QgsRasterCalculator::Result
Result
Result of the calculation.
Definition: qgsrastercalculator.h:81
QgsRasterLayer::bandCount
int bandCount() const
Returns the number of bands in this layer.
Definition: qgsrasterlayer.cpp:216
outputCrs
const QgsCoordinateReferenceSystem & outputCrs
Definition: qgswfsgetfeature.cpp:115
QgsCoordinateTransformContext
Definition: qgscoordinatetransformcontext.h:57
qgsrasterlayer.h
QgsRasterCalculatorEntry::ref
QString ref
Name of entry.
Definition: qgsrastercalculator.h:59
QgsRasterMatrix::isNumber
bool isNumber() const
Returns true if matrix is 1x1 (=scalar number)
Definition: qgsrastermatrix.h:75
QgsRasterMatrix::number
double number() const
Definition: qgsrastermatrix.h:76
QgsRasterCalcNode
Definition: qgsrastercalcnode.h:35
QgsCoordinateReferenceSystem::WKT_PREFERRED_GDAL
@ WKT_PREFERRED_GDAL
Preferred format for conversion of CRS to WKT for use with the GDAL library.
Definition: qgscoordinatereferencesystem.h:680
QgsProject::mapLayers
QMap< QString, QgsMapLayer * > mapLayers(const bool validOnly=false) const
Returns a map of all registered layers by layer ID.
Definition: qgsproject.cpp:3347
QgsFeedback::canceled
void canceled()
Internal routines can connect to this signal if they use event loop.
qgsgdalutils.h
QgsRasterCalculator::Canceled
@ Canceled
User canceled calculation.
Definition: qgsrastercalculator.h:86
QgsProject::transformContext
QgsCoordinateTransformContext transformContext
Definition: qgsproject.h:99
QgsRasterCalculator::BandError
@ BandError
Invalid band number for input.
Definition: qgsrastercalculator.h:89
QgsRasterMatrix::data
double * data()
Returns data array (but not ownership)
Definition: qgsrastermatrix.h:82
QgsProject::instance
static QgsProject * instance()
Returns the QgsProject singleton instance.
Definition: qgsproject.cpp:458
qgsogrutils.h
QgsGdalUtils::supportsRasterCreate
static bool supportsRasterCreate(GDALDriverH driver)
Reads whether a driver supports GDALCreate() for raster purposes.
Definition: qgsgdalutils.cpp:30
QgsFeedback::cancel
void cancel()
Tells the internal routines that the current operation should be canceled. This should be run by the ...
Definition: qgsfeedback.h:97
QgsDebugMsg
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
QgsRasterCalculatorEntry
Definition: qgsrastercalculator.h:39
QgsRectangle
Definition: qgsrectangle.h:41
QgsOpenClUtils::available
static bool available()
Checks whether a suitable OpenCL platform and device is available on this system and initialize the Q...
Definition: qgsopenclutils.cpp:463
QgsRasterMatrix::setNodataValue
void setNodataValue(double d)
Definition: qgsrastermatrix.h:96
QgsRasterProjector
QgsRasterProjector implements approximate projection support for it calculates grid of points in sour...
Definition: qgsrasterprojector.h:47
QgsRasterDataProvider::block
QgsRasterBlock * block(int bandNo, const QgsRectangle &boundingBox, int width, int height, QgsRasterBlockFeedback *feedback=nullptr) override
Read block of data using given extent and size.
Definition: qgsrasterdataprovider.cpp:45
qgsrasterinterface.h
QgsOpenClUtils::CPLAllocator
Tiny smart-pointer-like wrapper around CPLMalloc and CPLFree: this is needed because OpenCL C++ API m...
Definition: qgsopenclutils.h:237
QgsRasterCalculator::processCalculation
Result processCalculation(QgsFeedback *feedback=nullptr)
Starts the calculation and writes a new raster.
Definition: qgsrastercalculator.cpp:100
qgsrastermatrix.h
QgsRasterInterface::setInput
virtual bool setInput(QgsRasterInterface *input)
Set input.
Definition: qgsrasterinterface.h:266
gdal::dataset_unique_ptr
std::unique_ptr< std::remove_pointer< GDALDatasetH >::type, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
Definition: qgsogrutils.h:134
QgsFeedback
Definition: qgsfeedback.h:43
QgsRasterCalculator::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)
QgsRasterCalculator constructor.
Definition: qgsrastercalculator.cpp:39
QgsCoordinateReferenceSystem::toWkt
QString toWkt(WktVariant variant=WKT1_GDAL, bool multiline=false, int indentationWidth=4) const
Returns a WKT representation of this CRS.
Definition: qgscoordinatereferencesystem.cpp:1931
QgsRasterProjector::setCrs
Q_DECL_DEPRECATED void setCrs(const QgsCoordinateReferenceSystem &srcCRS, const QgsCoordinateReferenceSystem &destCRS, int srcDatumTransform=-1, int destDatumTransform=-1)
Sets the source and destination CRS.
QgsOpenClUtils::enabled
static bool enabled()
Returns true if OpenCL is enabled in the user settings.
Definition: qgsopenclutils.cpp:265
typeName
const QString & typeName
Definition: qgswfsgetfeature.cpp:109
gdal::fast_delete_and_close
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:71
QgsRasterLayer
Definition: qgsrasterlayer.h:72
QgsRasterCalculator::CreateOutputError
@ CreateOutputError
Error creating output data file.
Definition: qgsrastercalculator.h:84
QgsCoordinateReferenceSystem
Definition: qgscoordinatereferencesystem.h:206
QgsRasterCalculator::CalculationError
@ CalculationError
Error occurred while performing calculation.
Definition: qgsrastercalculator.h:90
QgsRectangle::yMaximum
double yMaximum() const
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:172
QgsRasterProjector::setPrecision
void setPrecision(Precision precision)
Definition: qgsrasterprojector.h:93
QgsRasterCalculator::lastError
QString lastError() const
Returns a description of the last error encountered.
Definition: qgsrastercalculator.cpp:677
QgsFeedback::isCanceled
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:66
QgsRasterMatrix
Definition: qgsrastermatrix.h:28
QgsRasterCalculator::Success
@ Success
Calculation successful.
Definition: qgsrastercalculator.h:83
QgsRasterCalculator::InputLayerError
@ InputLayerError
Error reading input layer.
Definition: qgsrastercalculator.h:85
QgsRasterLayer::providerType
QString providerType() const
[ data provider interface ] Which provider is being used for this Raster Layer?
Definition: qgsrasterlayer.cpp:545
QgsRasterCalculator::MemoryError
@ MemoryError
Error allocating memory for result.
Definition: qgsrastercalculator.h:88
QgsRasterCalculator::ParserError
@ ParserError
Error parsing formula.
Definition: qgsrastercalculator.h:87
QgsRasterProjector::Exact
@ Exact
Exact, precise but slow.
Definition: qgsrasterprojector.h:60
QgsRasterBlockFeedback
Definition: qgsrasterinterface.h:40
QgsRectangle::height
double height() const
Returns the height of the rectangle.
Definition: qgsrectangle.h:209
QgsRasterCalculatorEntry::rasterEntries
static QVector< QgsRasterCalculatorEntry > rasterEntries()
Creates a list of raster entries from the current project.
Definition: qgsrastercalculator.cpp:682
QgsOpenClUtils::commandQueue
static cl::CommandQueue commandQueue()
Create an OpenCL command queue from the default context.
Definition: qgsopenclutils.cpp:607
qgsrastercalculator.h
QgsOpenClUtils::context
static cl::Context context()
Context factory.
Definition: qgsopenclutils.cpp:628
QgsRasterCalcNode::parseRasterCalcString
static QgsRasterCalcNode * parseRasterCalcString(const QString &str, QString &parserErrorMsg)
Definition: qgsrastercalcnode.cpp:379
qgsopenclutils.h
QgsRasterProjector::block
QgsRasterBlock * block(int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback=nullptr) override
Read block of data using given extent and size.
Definition: qgsrasterprojector.cpp:771
qgsfeedback.h
QgsOpenClUtils::buildProgram
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...
Definition: qgsopenclutils.cpp:642
qgsproject.h
QgsRectangle::width
double width() const
Returns the width of the rectangle.
Definition: qgsrectangle.h:202
QgsRectangle::setYMinimum
void setYMinimum(double y)
Set the minimum y value.
Definition: qgsrectangle.h:140
QgsRectangle::setYMaximum
void setYMaximum(double y)
Set the maximum y value.
Definition: qgsrectangle.h:145
QgsRasterLayer::dataProvider
QgsRasterDataProvider * dataProvider() override
Returns the source data provider.
Definition: qgsrasterlayer.cpp:233
QgsRasterCalculatorEntry::bandNumber
int bandNumber
Band number for entry.
Definition: qgsrastercalculator.h:69
QgsRectangle::xMinimum
double xMinimum() const
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:167
qgsrasterdataprovider.h