QGIS API Documentation  3.8.0-Zanzibar (11aff65)
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 #ifdef HAVE_OPENCL
128  // Check for matrix nodes, GPU implementation does not support them
129  QList<const QgsRasterCalcNode *> nodeList;
130  if ( QgsOpenClUtils::enabled() && QgsOpenClUtils::available() && calcNode->findNodes( QgsRasterCalcNode::Type::tMatrix ).isEmpty() )
131  return processCalculationGPU( std::move( calcNode ), feedback );
132 #endif
133 
134  //open output dataset for writing
135  GDALDriverH outputDriver = openOutputDriver();
136  if ( !outputDriver )
137  {
138  mLastError = QObject::tr( "Could not obtain driver for %1" ).arg( mOutputFormat );
139  return CreateOutputError;
140  }
141 
142  gdal::dataset_unique_ptr outputDataset( openOutputFile( outputDriver ) );
143  if ( !outputDataset )
144  {
145  mLastError = QObject::tr( "Could not create output %1" ).arg( mOutputFile );
146  return CreateOutputError;
147  }
148 
149  GDALSetProjection( outputDataset.get(), mOutputCrs.toWkt().toLocal8Bit().data() );
150  GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
151 
152  float outputNodataValue = -FLT_MAX;
153  GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
154 
155  // Check if we need to read the raster as a whole (which is memory inefficient
156  // and not interruptable by the user) by checking if any raster matrix nodes are
157  // in the expression
158  bool requiresMatrix = ! calcNode->findNodes( QgsRasterCalcNode::Type::tMatrix ).isEmpty();
159 
160  // Take the fast route (process one line at a time) if we can
161  if ( ! requiresMatrix )
162  {
163  // Map of raster names -> blocks
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 ) )
167  {
168  QString layerRef( r->toString().remove( 0, 1 ) );
169  layerRef.chop( 1 );
170  if ( ! inputBlocks.count( layerRef ) )
171  {
172  for ( const auto &ref : mRasterEntries )
173  {
174  if ( ref.ref == layerRef )
175  {
176  uniqueRasterEntries[layerRef] = ref;
177  inputBlocks[layerRef ] = qgis::make_unique<QgsRasterBlock>();
178  }
179  }
180  }
181  }
182 
183  //read / write line by line
184  QMap<QString, QgsRasterBlock * > _rasterData;
185  // Cast to float
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 )
189  {
190  if ( feedback )
191  {
192  feedback->setProgress( 100.0 * static_cast< double >( row ) / mNumOutputRows );
193  }
194 
195  if ( feedback && feedback->isCanceled() )
196  {
197  break;
198  }
199 
200  // Calculates the rect for a single row read
201  QgsRectangle rect( mOutputRectangle );
202  rect.setYMaximum( rect.yMaximum() - rowHeight * row );
203  rect.setYMinimum( rect.yMaximum() - rowHeight );
204 
205  // Read rows into input blocks
206  for ( auto &layerRef : inputBlocks )
207  {
208  QgsRasterCalculatorEntry ref = uniqueRasterEntries[layerRef.first];
209  if ( uniqueRasterEntries[layerRef.first].raster->crs() != mOutputCrs )
210  {
211  QgsRasterProjector proj;
212  proj.setCrs( ref.raster->crs(), mOutputCrs, mTransformContext );
213  proj.setInput( ref.raster->dataProvider() );
215  layerRef.second.reset( proj.block( ref.bandNumber, rect, mNumOutputColumns, 1 ) );
216  }
217  else
218  {
219  inputBlocks[layerRef.first].reset( ref.raster->dataProvider()->block( ref.bandNumber, rect, mNumOutputColumns, 1 ) );
220  }
221  }
222 
223  QgsRasterMatrix resultMatrix;
224  resultMatrix.setNodataValue( outputNodataValue );
225 
226  _rasterData.clear();
227  for ( const auto &layerRef : inputBlocks )
228  {
229  _rasterData.insert( layerRef.first, inputBlocks[layerRef.first].get() );
230  }
231 
232  if ( calcNode->calculate( _rasterData, resultMatrix, 0 ) )
233  {
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 )
236  {
237  QgsDebugMsg( QStringLiteral( "RasterIO error!" ) );
238  }
239  }
240  }
241 
242  if ( feedback )
243  {
244  feedback->setProgress( 100.0 );
245  }
246  }
247  else // Original code (memory inefficient route)
248  {
249  QMap< QString, QgsRasterBlock * > inputBlocks;
250  QVector<QgsRasterCalculatorEntry>::const_iterator it = mRasterEntries.constBegin();
251  for ( ; it != mRasterEntries.constEnd(); ++it )
252  {
253 
254  std::unique_ptr< QgsRasterBlock > block;
255  // if crs transform needed
256  if ( it->raster->crs() != mOutputCrs )
257  {
258  QgsRasterProjector proj;
259  proj.setCrs( it->raster->crs(), mOutputCrs, it->raster->transformContext() );
260  proj.setInput( it->raster->dataProvider() );
262 
263  QgsRasterBlockFeedback *rasterBlockFeedback = new QgsRasterBlockFeedback();
264  QObject::connect( feedback, &QgsFeedback::canceled, rasterBlockFeedback, &QgsRasterBlockFeedback::cancel );
265  block.reset( proj.block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows, rasterBlockFeedback ) );
266  if ( rasterBlockFeedback->isCanceled() )
267  {
268  qDeleteAll( inputBlocks );
269  return Canceled;
270  }
271  }
272  else
273  {
274  block.reset( it->raster->dataProvider()->block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows ) );
275  }
276  if ( block->isEmpty() )
277  {
278  mLastError = QObject::tr( "Could not allocate required memory for %1" ).arg( it->ref );
279  qDeleteAll( inputBlocks );
280  return MemoryError;
281  }
282  inputBlocks.insert( it->ref, block.release() );
283  }
284 
285  QgsRasterMatrix resultMatrix;
286  resultMatrix.setNodataValue( outputNodataValue );
287 
288  //read / write line by line
289  for ( int i = 0; i < mNumOutputRows; ++i )
290  {
291  if ( feedback )
292  {
293  feedback->setProgress( 100.0 * static_cast< double >( i ) / mNumOutputRows );
294  }
295 
296  if ( feedback && feedback->isCanceled() )
297  {
298  break;
299  }
300 
301  if ( calcNode->calculate( inputBlocks, resultMatrix, i ) )
302  {
303  bool resultIsNumber = resultMatrix.isNumber();
304  float *calcData = new float[mNumOutputColumns];
305 
306  for ( int j = 0; j < mNumOutputColumns; ++j )
307  {
308  calcData[j] = ( float )( resultIsNumber ? resultMatrix.number() : resultMatrix.data()[j] );
309  }
310 
311  //write scanline to the dataset
312  if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
313  {
314  QgsDebugMsg( QStringLiteral( "RasterIO error!" ) );
315  }
316 
317  delete[] calcData;
318  }
319 
320  }
321 
322  if ( feedback )
323  {
324  feedback->setProgress( 100.0 );
325  }
326 
327  //close datasets and release memory
328  calcNode.reset();
329  qDeleteAll( inputBlocks );
330  inputBlocks.clear();
331 
332  }
333 
334  if ( feedback && feedback->isCanceled() )
335  {
336  //delete the dataset without closing (because it is faster)
337  gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
338  return Canceled;
339  }
340  return Success;
341 }
342 
343 #ifdef HAVE_OPENCL
344 QgsRasterCalculator::Result QgsRasterCalculator::processCalculationGPU( std::unique_ptr< QgsRasterCalcNode > calcNode, QgsFeedback *feedback )
345 {
346 
347  QString cExpression( calcNode->toString( true ) );
348 
349  QList<const QgsRasterCalcNode *> nodeList( calcNode->findNodes( QgsRasterCalcNode::Type::tRasterRef ) );
350  QSet<QString> capturedTexts;
351  for ( const auto &r : qgis::as_const( nodeList ) )
352  {
353  QString s( r->toString().remove( 0, 1 ) );
354  s.chop( 1 );
355  capturedTexts.insert( s );
356  }
357 
358  // Extract all references
359  struct LayerRef
360  {
361  QString name;
362  int band;
363  QgsRasterLayer *layer = nullptr;
364  QString varName;
365  QString typeName;
366  size_t index;
367  size_t bufferSize;
368  size_t dataSize;
369  };
370 
371  // Collects all layers, band, name, varName and size information
372  std::vector<LayerRef> inputRefs;
373  size_t refCounter = 0;
374  for ( const auto &r : capturedTexts )
375  {
376  if ( r.startsWith( '"' ) )
377  continue;
378  QStringList parts( r.split( '@' ) );
379  if ( parts.count() != 2 )
380  continue;
381  bool ok = false;
382  LayerRef entry;
383  entry.name = r;
384  entry.band = parts[1].toInt( &ok );
385  for ( const auto &ref : mRasterEntries )
386  {
387  if ( ref.ref == entry.name )
388  entry.layer = ref.raster;
389  }
390  if ( !( entry.layer && entry.layer->dataProvider() && ok ) )
391  continue;
392  entry.dataSize = entry.layer->dataProvider()->dataTypeSize( entry.band );
393  switch ( entry.layer->dataProvider()->dataType( entry.band ) )
394  {
395  case Qgis::DataType::Byte:
396  entry.typeName = QStringLiteral( "unsigned char" );
397  break;
398  case Qgis::DataType::UInt16:
399  entry.typeName = QStringLiteral( "unsigned int" );
400  break;
401  case Qgis::DataType::Int16:
402  entry.typeName = QStringLiteral( "short" );
403  break;
404  case Qgis::DataType::UInt32:
405  entry.typeName = QStringLiteral( "unsigned int" );
406  break;
407  case Qgis::DataType::Int32:
408  entry.typeName = QStringLiteral( "int" );
409  break;
410  case Qgis::DataType::Float32:
411  entry.typeName = QStringLiteral( "float" );
412  break;
413  // FIXME: not sure all OpenCL implementations support double
414  // maybe safer to fall back to the CPU implementation
415  // after a compatibility check
416  case Qgis::DataType::Float64:
417  entry.typeName = QStringLiteral( "double" );
418  break;
419  default:
420  return BandError;
421  }
422  entry.bufferSize = entry.dataSize * mNumOutputColumns;
423  entry.index = refCounter;
424  entry.varName = QStringLiteral( "input_raster_%1_band_%2" )
425  .arg( refCounter++ )
426  .arg( entry.band );
427  inputRefs.push_back( entry );
428  }
429 
430  // Prepare context and queue
431  cl::Context ctx( QgsOpenClUtils::context() );
432  cl::CommandQueue queue( QgsOpenClUtils::commandQueue() );
433 
434  // Create the C expression
435  std::vector<cl::Buffer> inputBuffers;
436  inputBuffers.reserve( inputRefs.size() );
437  QStringList inputArgs;
438  for ( const auto &ref : inputRefs )
439  {
440  cExpression.replace( QStringLiteral( "\"%1\"" ).arg( ref.name ), QStringLiteral( "%1[i]" ).arg( ref.varName ) );
441  inputArgs.append( QStringLiteral( "__global %1 *%2" )
442  .arg( ref.typeName )
443  .arg( ref.varName ) );
444  inputBuffers.push_back( cl::Buffer( ctx, CL_MEM_READ_ONLY, ref.bufferSize, nullptr, nullptr ) );
445  }
446 
447  //qDebug() << cExpression;
448 
449  // Create the program
450  QString programTemplate( R"CL(
451  // Inputs:
452  ##INPUT_DESC##
453  // Expression: ##EXPRESSION_ORIGINAL##
454  __kernel void rasterCalculator( ##INPUT##
455  __global float *resultLine
456  )
457  {
458  // Get the index of the current element
459  const int i = get_global_id(0);
460  // Expression
461  resultLine[i] = ##EXPRESSION##;
462  }
463  )CL" );
464 
465  QStringList inputDesc;
466  for ( const auto &ref : inputRefs )
467  {
468  inputDesc.append( QStringLiteral( " // %1 = %2" ).arg( ref.varName ).arg( ref.name ) );
469  }
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( ) );
474 
475  // qDebug() << programTemplate;
476 
477  // Create a program from the kernel source
478  cl::Program program( QgsOpenClUtils::buildProgram( programTemplate, QgsOpenClUtils::ExceptionBehavior::Throw ) );
479 
480  // Create the buffers, output is float32 (4 bytes)
481  // We assume size of float = 4 because that's the size used by OpenCL and IEEE 754
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 );
486 
487  auto kernel = cl::Kernel( program, "rasterCalculator" );
488 
489  for ( unsigned int i = 0; i < inputBuffers.size() ; i++ )
490  {
491  kernel.setArg( i, inputBuffers.at( i ) );
492  }
493  kernel.setArg( static_cast<unsigned int>( inputBuffers.size() ), resultLineBuffer );
494 
495  QgsOpenClUtils::CPLAllocator<float> resultLine( static_cast<size_t>( mNumOutputColumns ) );
496 
497  //open output dataset for writing
498  GDALDriverH outputDriver = openOutputDriver();
499  if ( !outputDriver )
500  {
501  mLastError = QObject::tr( "Could not obtain driver for %1" ).arg( mOutputFormat );
502  return CreateOutputError;
503  }
504 
505  gdal::dataset_unique_ptr outputDataset( openOutputFile( outputDriver ) );
506  if ( !outputDataset )
507  {
508  mLastError = QObject::tr( "Could not create output %1" ).arg( mOutputFile );
509  return CreateOutputError;
510  }
511 
512  GDALSetProjection( outputDataset.get(), mOutputCrs.toWkt().toLocal8Bit().data() );
513 
514  GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
515  if ( !outputRasterBand )
516  return BandError;
517 
518  // Input block (buffer)
519  std::unique_ptr<QgsRasterBlock> block;
520 
521  // Run kernel on all scanlines
522  auto rowHeight = mOutputRectangle.height() / mNumOutputRows;
523  for ( int line = 0; line < mNumOutputRows; line++ )
524  {
525  if ( feedback && feedback->isCanceled() )
526  {
527  break;
528  }
529 
530  if ( feedback )
531  {
532  feedback->setProgress( 100.0 * static_cast< double >( line ) / mNumOutputRows );
533  }
534 
535  // Read lines from rasters into the buffers
536  for ( const auto &ref : inputRefs )
537  {
538  // Read one row
539  QgsRectangle rect( mOutputRectangle );
540  rect.setYMaximum( rect.yMaximum() - rowHeight * line );
541  rect.setYMinimum( rect.yMaximum() - rowHeight );
542 
543  // TODO: check if this is too slow
544  // if crs transform needed
545  if ( ref.layer->crs() != mOutputCrs )
546  {
547  QgsRasterProjector proj;
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 ) );
552  }
553  else
554  {
555  block.reset( ref.layer->dataProvider()->block( ref.band, rect, mNumOutputColumns, 1 ) );
556  }
557 
558  //for ( int i = 0; i < mNumOutputColumns; i++ )
559  // qDebug() << "Input: " << line << i << ref.varName << " = " << block->value( 0, i );
560  //qDebug() << "Writing buffer " << ref.index;
561 
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() );
565 
566  }
567  // Run the kernel
568  queue.enqueueNDRangeKernel(
569  kernel,
570  0,
571  cl::NDRange( mNumOutputColumns )
572  );
573 
574  // Write the result
575  queue.enqueueReadBuffer( resultLineBuffer, CL_TRUE, 0,
576  resultBufferSize, resultLine.get() );
577 
578  //for ( int i = 0; i < mNumOutputColumns; i++ )
579  // qDebug() << "Output: " << line << i << " = " << resultLine[i];
580 
581  if ( GDALRasterIO( outputRasterBand, GF_Write, 0, line, mNumOutputColumns, 1, resultLine.get(), mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
582  {
583  return CreateOutputError;
584  }
585  }
586 
587  if ( feedback && feedback->isCanceled() )
588  {
589  //delete the dataset without closing (because it is faster)
590  gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
591  return Canceled;
592  }
593 
594  inputBuffers.clear();
595 
596  return Success;
597 }
598 #endif
599 
600 GDALDriverH QgsRasterCalculator::openOutputDriver()
601 {
602  //open driver
603  GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
604 
605  if ( !outputDriver )
606  {
607  return outputDriver; //return nullptr, driver does not exist
608  }
609 
610  if ( !QgsGdalUtils::supportsRasterCreate( outputDriver ) )
611  {
612  return nullptr; //driver exist, but it does not support the create operation
613  }
614 
615  return outputDriver;
616 }
617 
618 gdal::dataset_unique_ptr QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver )
619 {
620  //open output file
621  char **papszOptions = nullptr;
622  gdal::dataset_unique_ptr outputDataset( GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions ) );
623  if ( !outputDataset )
624  {
625  return nullptr;
626  }
627 
628  //assign georef information
629  double geotransform[6];
630  outputGeoTransform( geotransform );
631  GDALSetGeoTransform( outputDataset.get(), geotransform );
632 
633  return outputDataset;
634 }
635 
636 void QgsRasterCalculator::outputGeoTransform( double *transform ) const
637 {
638  transform[0] = mOutputRectangle.xMinimum();
639  transform[1] = mOutputRectangle.width() / mNumOutputColumns;
640  transform[2] = 0;
641  transform[3] = mOutputRectangle.yMaximum();
642  transform[4] = 0;
643  transform[5] = -mOutputRectangle.height() / mNumOutputRows;
644 }
645 
647 {
648  return mLastError;
649 }
650 
651 QVector<QgsRasterCalculatorEntry> QgsRasterCalculatorEntry::rasterEntries()
652 {
653  QVector<QgsRasterCalculatorEntry> availableEntries;
654  const QMap<QString, QgsMapLayer *> &layers = QgsProject::instance()->mapLayers();
655 
656  auto uniqueRasterBandIdentifier = [ & ]( QgsRasterCalculatorEntry & entry ) -> bool
657  {
658  unsigned int i( 1 );
659  entry.ref = QStringLiteral( "%1@%2" ).arg( entry.raster->name() ).arg( entry.bandNumber );
660  while ( true )
661  {
662  bool unique( true );
663  for ( const auto &ref : qgis::as_const( availableEntries ) )
664  {
665  // Safety belt
666  if ( !( entry.raster && ref.raster ) )
667  continue;
668  // Check if is another band of the same raster
669  if ( ref.raster->publicSource() == entry.raster->publicSource() )
670  {
671  if ( ref.bandNumber != entry.bandNumber )
672  {
673  continue;
674  }
675  else // a layer with the same data source was already added to the list
676  {
677  return false;
678  }
679  }
680  // If same name but different source
681  if ( ref.ref == entry.ref )
682  {
683  unique = false;
684  entry.ref = QStringLiteral( "%1_%2@%3" ).arg( entry.raster->name() ).arg( i++ ).arg( entry.bandNumber );
685  }
686  }
687  if ( unique )
688  return true;
689  }
690  };
691 
692  QMap<QString, QgsMapLayer *>::const_iterator layerIt = layers.constBegin();
693  for ( ; layerIt != layers.constEnd(); ++layerIt )
694  {
695  QgsRasterLayer *rlayer = qobject_cast<QgsRasterLayer *>( layerIt.value() );
696  if ( rlayer && rlayer->dataProvider() && rlayer->dataProvider()->name() == QLatin1String( "gdal" ) )
697  {
698  //get number of bands
699  for ( int i = 0; i < rlayer->bandCount(); ++i )
700  {
702  entry.raster = rlayer;
703  entry.bandNumber = i + 1;
704  if ( ! uniqueRasterBandIdentifier( entry ) )
705  break;
706  availableEntries.push_back( entry );
707  }
708  }
709  }
710  return availableEntries;
711 }
A rectangle specified with double values.
Definition: qgsrectangle.h:41
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 ...
Definition: qgsfeedback.h:85
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
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.
Definition: qgsfeedback.h:63
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...
Definition: qgsfeedback.h:44
QgsRasterDataProvider * dataProvider() override
Returns the layer&#39;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...
Definition: qgsogrutils.cpp:62
Calculation successful.
const QString & typeName
void setYMinimum(double y)
Set the minimum y value.
Definition: qgsrectangle.h:140
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
Definition: qgsproject.h:96
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.
Definition: qgsfeedback.h:54
double number() const
void setYMaximum(double y)
Set the maximum y value.
Definition: qgsrectangle.h:145
static QgsProject * instance()
Returns the QgsProject singleton instance.
Definition: qgsproject.cpp:438
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).
Definition: qgsrectangle.h:172
Exact, precise but slow.
std::unique_ptr< std::remove_pointer< GDALDatasetH >::type, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
Definition: qgsogrutils.h:134
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
Definition: qgsmaplayer.h:85
double height() const
Returns the height of the rectangle.
Definition: qgsrectangle.h:209
Q_DECL_DEPRECATED void setCrs(const QgsCoordinateReferenceSystem &srcCRS, const QgsCoordinateReferenceSystem &destCRS, int srcDatumTransform=-1, int destDatumTransform=-1)
Sets the source and destination CRS.