Statistics
| Branch: | Tag: | Revision:

qgis / src / analysis / raster / qgsrastercalculator.cpp @ master

History | View | Annotate | Download (26.5 KB)

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

    
135
QgsRasterCalculator::Result QgsRasterCalculator::processCalculation( QgsFeedback *feedback )
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() );
252
          proj.setPrecision( QgsRasterProjector::Exact );
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() );
305
        proj.setPrecision( QgsRasterProjector::Exact );
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
    {
449
      case Qgis::DataType::Byte:
450
        entry.typeName = QStringLiteral( "unsigned char" );
451
        break;
452
      case Qgis::DataType::UInt16:
453
        entry.typeName = QStringLiteral( "unsigned int" );
454
        break;
455
      case Qgis::DataType::Int16:
456
        entry.typeName = QStringLiteral( "short" );
457
        break;
458
      case Qgis::DataType::UInt32:
459
        entry.typeName = QStringLiteral( "unsigned int" );
460
        break;
461
      case Qgis::DataType::Int32:
462
        entry.typeName = QStringLiteral( "int" );
463
        break;
464
      case Qgis::DataType::Float32:
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
470
      case Qgis::DataType::Float64:
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() );
621
          proj.setPrecision( QgsRasterProjector::Exact );
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

    
726
QString QgsRasterCalculator::lastError() const
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
      {
781
        QgsRasterCalculatorEntry entry;
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
}