QGIS API Documentation 3.40.0-Bratislava (b56115d8743)
Loading...
Searching...
No Matches
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"
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//
43int 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
74QgsRasterCalculator::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
87QgsRasterCalculator::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!
104QgsRasterCalculator::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!
121QgsRasterCalculator::QgsRasterCalculator( const QString &formulaString, const QString &outputFile, const QString &outputFormat,
122 const QgsRectangle &outputExtent, const QgsCoordinateReferenceSystem &outputCrs, int nOutputColumns, int nOutputRows, const QVector<QgsRasterCalculatorEntry> &rasterEntries )
123 : mFormulaString( formulaString )
124 , mOutputFile( outputFile )
125 , mOutputFormat( outputFormat )
126 , mOutputRectangle( outputExtent )
127 , mOutputCrs( outputCrs )
128 , mNumOutputColumns( nOutputColumns )
129 , mNumOutputRows( nOutputRows )
130 , mRasterEntries( rasterEntries )
131 , mTransformContext( QgsProject::instance()->transformContext() )
132{
133}
134
136{
137 mLastError.clear();
138
139 //prepare search string / tree
140 std::unique_ptr< QgsRasterCalcNode > calcNode( QgsRasterCalcNode::parseRasterCalcString( mFormulaString, mLastError ) );
141 if ( !calcNode )
142 {
143 //error
144 return ParserError;
145 }
146
147 // Check input layers and bands
148 for ( const auto &entry : std::as_const( mRasterEntries ) )
149 {
150 if ( !entry.raster ) // no raster layer in entry
151 {
152 mLastError = QObject::tr( "No raster layer for entry %1" ).arg( entry.ref );
153 return InputLayerError;
154 }
155 if ( entry.bandNumber <= 0 || entry.bandNumber > entry.raster->bandCount() )
156 {
157 mLastError = QObject::tr( "Band number %1 is not valid for entry %2" ).arg( entry.bandNumber ).arg( entry.ref );
158 return BandError;
159 }
160 }
161
162 // Check if we need to read the raster as a whole (which is memory inefficient
163 // and not interruptible by the user) by checking if any raster matrix nodes are
164 // in the expression
165 bool requiresMatrix = ! calcNode->findNodes( QgsRasterCalcNode::Type::tMatrix ).isEmpty();
166
167#ifdef HAVE_OPENCL
168 // Check for matrix nodes, GPU implementation does not support them
169 if ( QgsOpenClUtils::enabled() && QgsOpenClUtils::available() && ! requiresMatrix )
170 {
171 return processCalculationGPU( std::move( calcNode ), feedback );
172 }
173#endif
174
175 //open output dataset for writing
176 GDALDriverH outputDriver = openOutputDriver();
177 if ( !outputDriver )
178 {
179 mLastError = QObject::tr( "Could not obtain driver for %1" ).arg( mOutputFormat );
180 return CreateOutputError;
181 }
182
183 gdal::dataset_unique_ptr outputDataset( openOutputFile( outputDriver ) );
184 if ( !outputDataset )
185 {
186 mLastError = QObject::tr( "Could not create output %1" ).arg( mOutputFile );
187 return CreateOutputError;
188 }
189
190 GDALSetProjection( outputDataset.get(), mOutputCrs.toWkt( Qgis::CrsWktVariant::PreferredGdal ).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 {
250 proj.setCrs( ref.raster->crs(), mOutputCrs, mTransformContext );
251 proj.setInput( ref.raster->dataProvider() );
253 layerRef.second.reset( proj.block( ref.bandNumber, rect, mNumOutputColumns, 1 ) );
254 }
255 else
256 {
257 layerRef.second.reset( ref.raster->dataProvider()->block( ref.bandNumber, rect, mNumOutputColumns, 1 ) );
258 }
259 }
260
261 // 1 row X mNumOutputColumns matrix
262 QgsRasterMatrix resultMatrix( mNumOutputColumns, 1, nullptr, outputNodataValue );
263
264 _rasterData.clear();
265 for ( const auto &layerRef : inputBlocks )
266 {
267 _rasterData.insert( layerRef.first, layerRef.second.get() );
268 }
269
270 if ( calcNode->calculate( _rasterData, resultMatrix, 0 ) )
271 {
272 std::copy( resultMatrix.data(), resultMatrix.data() + mNumOutputColumns, castedResult.begin() );
273 if ( GDALRasterIO( outputRasterBand, GF_Write, 0, row, mNumOutputColumns, 1, castedResult.data(), mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
274 {
275 QgsDebugError( 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 {
303 proj.setCrs( it->raster->crs(), mOutputCrs, it->raster->transformContext() );
304 proj.setInput( it->raster->dataProvider() );
306
307 QgsRasterBlockFeedback *rasterBlockFeedback = new QgsRasterBlockFeedback();
308 QObject::connect( feedback, &QgsFeedback::canceled, rasterBlockFeedback, &QgsRasterBlockFeedback::cancel );
309 block.reset( proj.block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows, rasterBlockFeedback ) );
310 if ( rasterBlockFeedback->isCanceled() )
311 {
312 qDeleteAll( inputBlocks );
313 return Canceled;
314 }
315 }
316 else
317 {
318 block.reset( it->raster->dataProvider()->block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows ) );
319 }
320 if ( block->isEmpty() )
321 {
322 mLastError = QObject::tr( "Could not allocate required memory for %1" ).arg( it->ref );
323 qDeleteAll( inputBlocks );
324 return MemoryError;
325 }
326 inputBlocks.insert( it->ref, block.release() );
327 }
328
329 QgsRasterMatrix resultMatrix;
330 resultMatrix.setNodataValue( outputNodataValue );
331
332 //read / write line by line
333 for ( int i = 0; i < mNumOutputRows; ++i )
334 {
335 if ( feedback )
336 {
337 feedback->setProgress( 100.0 * static_cast< double >( i ) / mNumOutputRows );
338 }
339
340 if ( feedback && feedback->isCanceled() )
341 {
342 break;
343 }
344
345 if ( calcNode->calculate( inputBlocks, resultMatrix, i ) )
346 {
347 bool resultIsNumber = resultMatrix.isNumber();
348 float *calcData = new float[mNumOutputColumns];
349
350 for ( int j = 0; j < mNumOutputColumns; ++j )
351 {
352 calcData[j] = ( float )( resultIsNumber ? resultMatrix.number() : resultMatrix.data()[j] );
353 }
354
355 //write scanline to the dataset
356 if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
357 {
358 QgsDebugError( 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
398QgsRasterCalculator::Result QgsRasterCalculator::processCalculationGPU( std::unique_ptr< QgsRasterCalcNode > calcNode, QgsFeedback *feedback )
399{
400
401 QString cExpression( calcNode->toString( true ) );
402
403 QList<const QgsRasterCalcNode *> nodeList( calcNode->findNodes( QgsRasterCalcNode::Type::tRasterRef ) );
404 QSet<QString> capturedTexts;
405 for ( const auto &r : std::as_const( nodeList ) )
406 {
407 QString s( r->toString().remove( 0, 1 ) );
408 s.chop( 1 );
409 capturedTexts.insert( s );
410 }
411
412 // Extract all references
413 struct LayerRef
414 {
415 QString name;
416 int band;
417 QgsRasterLayer *layer = nullptr;
418 QString varName;
419 QString typeName;
420 size_t index;
421 size_t bufferSize;
422 size_t dataSize;
423 };
424
425 // Collects all layers, band, name, varName and size information
426 std::vector<LayerRef> inputRefs;
427 size_t refCounter = 0;
428 for ( const auto &r : capturedTexts )
429 {
430 if ( r.startsWith( '"' ) )
431 continue;
432 QStringList parts( r.split( '@' ) );
433 if ( parts.count() != 2 )
434 continue;
435 bool ok = false;
436 LayerRef entry;
437 entry.name = r;
438 entry.band = parts[1].toInt( &ok );
439 for ( const auto &ref : std::as_const( mRasterEntries ) )
440 {
441 if ( ref.ref == entry.name )
442 entry.layer = ref.raster;
443 }
444 if ( !( entry.layer && entry.layer->dataProvider() && ok ) )
445 continue;
446 entry.dataSize = entry.layer->dataProvider()->dataTypeSize( entry.band );
447 switch ( entry.layer->dataProvider()->dataType( entry.band ) )
448 {
450 entry.typeName = QStringLiteral( "unsigned char" );
451 break;
453 entry.typeName = QStringLiteral( "signed char" );
454 break;
456 entry.typeName = QStringLiteral( "unsigned int" );
457 break;
459 entry.typeName = QStringLiteral( "short" );
460 break;
462 entry.typeName = QStringLiteral( "unsigned int" );
463 break;
465 entry.typeName = QStringLiteral( "int" );
466 break;
468 entry.typeName = QStringLiteral( "float" );
469 break;
470 // FIXME: not sure all OpenCL implementations support double
471 // maybe safer to fall back to the CPU implementation
472 // after a compatibility check
474 entry.typeName = QStringLiteral( "double" );
475 break;
483 return BandError;
484 }
485 entry.bufferSize = entry.dataSize * mNumOutputColumns;
486 entry.index = refCounter;
487 entry.varName = QStringLiteral( "input_raster_%1_band_%2" )
488 .arg( refCounter++ )
489 .arg( entry.band );
490 inputRefs.push_back( entry );
491 }
492
493 // May throw an openCL exception
494 try
495 {
496 // Prepare context and queue
497 cl::Context ctx( QgsOpenClUtils::context() );
498 cl::CommandQueue queue( QgsOpenClUtils::commandQueue() );
499
500 // Create the C expression
501 std::vector<cl::Buffer> inputBuffers;
502 inputBuffers.reserve( inputRefs.size() );
503 QStringList inputArgs;
504 for ( const auto &ref : inputRefs )
505 {
506 cExpression.replace( QStringLiteral( "\"%1\"" ).arg( ref.name ), QStringLiteral( "%1[i]" ).arg( ref.varName ) );
507 inputArgs.append( QStringLiteral( "__global %1 *%2" )
508 .arg( ref.typeName, ref.varName ) );
509 inputBuffers.push_back( cl::Buffer( ctx, CL_MEM_READ_ONLY, ref.bufferSize, nullptr, nullptr ) );
510 }
511
512 //qDebug() << cExpression;
513
514 // Create the program
515 QString programTemplate( R"CL(
516 // Inputs:
517 ##INPUT_DESC##
518 // Expression: ##EXPRESSION_ORIGINAL##
519 __kernel void rasterCalculator( ##INPUT##
520 __global float *resultLine
521 )
522 {
523 // Get the index of the current element
524 const int i = get_global_id(0);
525 // Check for nodata in input
526 if ( ##INPUT_NODATA_CHECK## )
527 resultLine[i] = -FLT_MAX;
528 // Expression
529 else
530 resultLine[i] = ##EXPRESSION##;
531 }
532 )CL" );
533
534 QStringList inputDesc;
535 QStringList inputNoDataCheck;
536 for ( const auto &ref : inputRefs )
537 {
538 inputDesc.append( QStringLiteral( " // %1 = %2" ).arg( ref.varName, ref.name ) );
539 if ( ref.layer->dataProvider()->sourceHasNoDataValue( ref.band ) )
540 {
541 inputNoDataCheck.append( QStringLiteral( "(float) %1[i] == (float) %2" ).arg( ref.varName, QString::number( ref.layer->dataProvider()->sourceNoDataValue( ref.band ), 'g', 10 ) ) );
542 }
543 }
544
545 programTemplate = programTemplate.replace( QLatin1String( "##INPUT_NODATA_CHECK##" ), inputNoDataCheck.isEmpty() ? QStringLiteral( "false" ) : inputNoDataCheck.join( QLatin1String( " || " ) ) );
546 programTemplate = programTemplate.replace( QLatin1String( "##INPUT_DESC##" ), inputDesc.join( '\n' ) );
547 programTemplate = programTemplate.replace( QLatin1String( "##INPUT##" ), !inputArgs.isEmpty() ? ( inputArgs.join( ',' ).append( ',' ) ) : QChar( ' ' ) );
548 programTemplate = programTemplate.replace( QLatin1String( "##EXPRESSION##" ), cExpression );
549 programTemplate = programTemplate.replace( QLatin1String( "##EXPRESSION_ORIGINAL##" ), calcNode->toString( ) );
550
551 //qDebug() << programTemplate;
552
553 // Create a program from the kernel source
554 cl::Program program( QgsOpenClUtils::buildProgram( programTemplate, QgsOpenClUtils::ExceptionBehavior::Throw ) );
555
556 // Create the buffers, output is float32 (4 bytes)
557 // We assume size of float = 4 because that's the size used by OpenCL and IEEE 754
558 Q_ASSERT( sizeof( float ) == 4 );
559 std::size_t resultBufferSize( 4 * static_cast<size_t>( mNumOutputColumns ) );
560 cl::Buffer resultLineBuffer( ctx, CL_MEM_WRITE_ONLY,
561 resultBufferSize, nullptr, nullptr );
562
563 auto kernel = cl::Kernel( program, "rasterCalculator" );
564
565 for ( unsigned int i = 0; i < inputBuffers.size() ; i++ )
566 {
567 kernel.setArg( i, inputBuffers.at( i ) );
568 }
569 kernel.setArg( static_cast<unsigned int>( inputBuffers.size() ), resultLineBuffer );
570
571 QgsOpenClUtils::CPLAllocator<float> resultLine( static_cast<size_t>( mNumOutputColumns ) );
572
573 //open output dataset for writing
574 GDALDriverH outputDriver = openOutputDriver();
575 if ( !outputDriver )
576 {
577 mLastError = QObject::tr( "Could not obtain driver for %1" ).arg( mOutputFormat );
578 return CreateOutputError;
579 }
580
581 gdal::dataset_unique_ptr outputDataset( openOutputFile( outputDriver ) );
582 if ( !outputDataset )
583 {
584 mLastError = QObject::tr( "Could not create output %1" ).arg( mOutputFile );
585 return CreateOutputError;
586 }
587
588 GDALSetProjection( outputDataset.get(), mOutputCrs.toWkt( Qgis::CrsWktVariant::PreferredGdal ).toLocal8Bit().data() );
589
590
591 GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
592 if ( !outputRasterBand )
593 return BandError;
594
595 const float outputNodataValue = -FLT_MAX;
596 GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );
597
598 // Input block (buffer)
599 std::unique_ptr<QgsRasterBlock> block;
600
601 // Run kernel on all scanlines
602 auto rowHeight = mOutputRectangle.height() / mNumOutputRows;
603 for ( int line = 0; line < mNumOutputRows; line++ )
604 {
605 if ( feedback && feedback->isCanceled() )
606 {
607 break;
608 }
609
610 if ( feedback )
611 {
612 feedback->setProgress( 100.0 * static_cast< double >( line ) / mNumOutputRows );
613 }
614
615 // Read lines from rasters into the buffers
616 for ( const auto &ref : inputRefs )
617 {
618 // Read one row
619 QgsRectangle rect( mOutputRectangle );
620 rect.setYMaximum( rect.yMaximum() - rowHeight * line );
621 rect.setYMinimum( rect.yMaximum() - rowHeight );
622
623 // TODO: check if this is too slow
624 // if crs transform needed
625 if ( ref.layer->crs() != mOutputCrs )
626 {
628 proj.setCrs( ref.layer->crs(), mOutputCrs, ref.layer->transformContext() );
629 proj.setInput( ref.layer->dataProvider() );
631 block.reset( proj.block( ref.band, rect, mNumOutputColumns, 1 ) );
632 }
633 else
634 {
635 block.reset( ref.layer->dataProvider()->block( ref.band, rect, mNumOutputColumns, 1 ) );
636 }
637
638 //for ( int i = 0; i < mNumOutputColumns; i++ )
639 // qDebug() << "Input: " << line << i << ref.varName << " = " << block->value( 0, i );
640 //qDebug() << "Writing buffer " << ref.index;
641
642 Q_ASSERT( ref.bufferSize == static_cast<size_t>( block->data().size( ) ) );
643 queue.enqueueWriteBuffer( inputBuffers[ref.index], CL_TRUE, 0,
644 ref.bufferSize, block->bits() );
645
646 }
647 // Run the kernel
648 queue.enqueueNDRangeKernel(
649 kernel,
650 0,
651 cl::NDRange( mNumOutputColumns )
652 );
653
654 // Write the result
655 queue.enqueueReadBuffer( resultLineBuffer, CL_TRUE, 0,
656 resultBufferSize, resultLine.get() );
657
658 //for ( int i = 0; i < mNumOutputColumns; i++ )
659 // qDebug() << "Output: " << line << i << " = " << resultLine[i];
660
661 if ( GDALRasterIO( outputRasterBand, GF_Write, 0, line, mNumOutputColumns, 1, resultLine.get(), mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
662 {
663 return CreateOutputError;
664 }
665 }
666
667 if ( feedback && feedback->isCanceled() )
668 {
669 //delete the dataset without closing (because it is faster)
670 gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
671 return Canceled;
672 }
673
674 inputBuffers.clear();
675
676 GDALComputeRasterStatistics( outputRasterBand, true, nullptr, nullptr, nullptr, nullptr, GdalProgressCallback, feedback );
677
678 }
679 catch ( cl::Error &e )
680 {
681 mLastError = e.what();
682 return CreateOutputError;
683 }
684
685 return Success;
686}
687#endif
688
689GDALDriverH QgsRasterCalculator::openOutputDriver()
690{
691 //open driver
692 GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
693
694 if ( !outputDriver )
695 {
696 return outputDriver; //return nullptr, driver does not exist
697 }
698
699 if ( !QgsGdalUtils::supportsRasterCreate( outputDriver ) )
700 {
701 return nullptr; //driver exist, but it does not support the create operation
702 }
703
704 return outputDriver;
705}
706
707gdal::dataset_unique_ptr QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver )
708{
709 //open output file
710 char **papszOptions = nullptr;
711 gdal::dataset_unique_ptr outputDataset( GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions ) );
712 if ( !outputDataset )
713 {
714 return nullptr;
715 }
716
717 //assign georef information
718 double geotransform[6];
719 outputGeoTransform( geotransform );
720 GDALSetGeoTransform( outputDataset.get(), geotransform );
721
722 return outputDataset;
723}
724
725void QgsRasterCalculator::outputGeoTransform( double *transform ) const
726{
727 transform[0] = mOutputRectangle.xMinimum();
728 transform[1] = mOutputRectangle.width() / mNumOutputColumns;
729 transform[2] = 0;
730 transform[3] = mOutputRectangle.yMaximum();
731 transform[4] = 0;
732 transform[5] = -mOutputRectangle.height() / mNumOutputRows;
733}
734
736{
737 return mLastError;
738}
739
740QVector<QgsRasterCalculatorEntry> QgsRasterCalculatorEntry::rasterEntries()
741{
742 QVector<QgsRasterCalculatorEntry> availableEntries;
743 const QMap<QString, QgsMapLayer *> &layers = QgsProject::instance()->mapLayers();
744
745 auto uniqueRasterBandIdentifier = [ & ]( QgsRasterCalculatorEntry & entry ) -> bool
746 {
747 unsigned int i( 1 );
748 entry.ref = QStringLiteral( "%1@%2" ).arg( entry.raster->name() ).arg( entry.bandNumber );
749 while ( true )
750 {
751 bool unique( true );
752 for ( const auto &ref : std::as_const( availableEntries ) )
753 {
754 // Safety belt
755 if ( !( entry.raster && ref.raster ) )
756 continue;
757 // Check if is another band of the same raster
758 if ( ref.raster->publicSource() == entry.raster->publicSource() )
759 {
760 if ( ref.bandNumber != entry.bandNumber )
761 {
762 continue;
763 }
764 else // a layer with the same data source was already added to the list
765 {
766 return false;
767 }
768 }
769 // If same name but different source
770 if ( ref.ref == entry.ref )
771 {
772 unique = false;
773 entry.ref = QStringLiteral( "%1_%2@%3" ).arg( entry.raster->name() ).arg( i++ ).arg( entry.bandNumber );
774 }
775 }
776 if ( unique )
777 return true;
778 }
779 };
780
781 QMap<QString, QgsMapLayer *>::const_iterator layerIt = layers.constBegin();
782 for ( ; layerIt != layers.constEnd(); ++layerIt )
783 {
784 QgsRasterLayer *rlayer = qobject_cast<QgsRasterLayer *>( layerIt.value() );
785 if ( rlayer && rlayer->dataProvider() && ( rlayer->dataProvider()->capabilities() & Qgis::RasterInterfaceCapability::Size ) )
786 {
787 //get number of bands
788 for ( int i = 0; i < rlayer->bandCount(); ++i )
789 {
791 entry.raster = rlayer;
792 entry.bandNumber = i + 1;
793 if ( ! uniqueRasterBandIdentifier( entry ) )
794 break;
795 availableEntries.push_back( entry );
796 }
797 }
798 }
799 return availableEntries;
800}
@ Size
Original data source size (and thus resolution) is known, it is not always available,...
@ CInt32
Complex Int32.
@ Float32
Thirty two bit floating point (float)
@ CFloat64
Complex Float64.
@ Int16
Sixteen bit signed integer (qint16)
@ ARGB32_Premultiplied
Color, alpha, red, green, blue, 4 bytes the same as QImage::Format_ARGB32_Premultiplied.
@ Int8
Eight bit signed integer (qint8) (added in QGIS 3.30)
@ UInt16
Sixteen bit unsigned integer (quint16)
@ Byte
Eight bit unsigned integer (quint8)
@ UnknownDataType
Unknown or unspecified type.
@ ARGB32
Color, alpha, red, green, blue, 4 bytes the same as QImage::Format_ARGB32.
@ Int32
Thirty two bit signed integer (qint32)
@ Float64
Sixty four bit floating point (double)
@ CFloat32
Complex Float32.
@ CInt16
Complex Int16.
@ UInt32
Thirty two bit unsigned integer (quint32)
@ PreferredGdal
Preferred format for conversion of CRS to WKT for use with the GDAL library.
This class represents a coordinate reference system (CRS).
QString toWkt(Qgis::CrsWktVariant variant=Qgis::CrsWktVariant::Wkt1Gdal, bool multiline=false, int indentationWidth=4) const
Returns a WKT representation of this CRS.
Contains information about the context in which a coordinate transform is executed.
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition qgsfeedback.h:44
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition qgsfeedback.h:53
void canceled()
Internal routines can connect to this signal if they use event loop.
void cancel()
Tells the internal routines that the current operation should be canceled. This should be run by the ...
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition qgsfeedback.h:61
static bool supportsRasterCreate(GDALDriverH driver)
Reads whether a driver supports GDALCreate() for raster purposes.
QgsCoordinateReferenceSystem crs
Definition qgsmaplayer.h:83
static Q_DECL_DEPRECATED cl::Program buildProgram(const cl::Context &context, const QString &source, ExceptionBehavior exceptionBehavior=Catch)
Build the program from source in the given context and depending on exceptionBehavior can throw or ca...
static cl::Context context()
Context factory.
static bool enabled()
Returns true if OpenCL is enabled in the user settings.
static bool available()
Checks whether a suitable OpenCL platform and device is available on this system and initialize the Q...
static cl::CommandQueue commandQueue()
Create an OpenCL command queue from the default context.
@ Throw
Write errors in the message log and re-throw exceptions.
Encapsulates a QGIS project, including sets of map layers and their styles, layouts,...
Definition qgsproject.h:107
static QgsProject * instance()
Returns the QgsProject singleton instance.
QgsCoordinateTransformContext transformContext
Definition qgsproject.h:113
QMap< QString, QgsMapLayer * > mapLayers(const bool validOnly=false) const
Returns a map of all registered layers by layer ID.
Feedback object tailored for raster block reading.
Represents a node in a raster calculator.
static QgsRasterCalcNode * parseRasterCalcString(const QString &str, QString &parserErrorMsg)
Represents an individual raster layer/band number entry within a raster calculation.
QgsRasterLayer * raster
Raster layer associated with entry.
static QVector< QgsRasterCalculatorEntry > rasterEntries()
Creates a list of raster entries from the current project.
int bandNumber
Band number for entry.
QString ref
Name of entry.
QgsRasterCalculator(const QString &formulaString, const QString &outputFile, const QString &outputFormat, const QgsRectangle &outputExtent, int nOutputColumns, int nOutputRows, const QVector< QgsRasterCalculatorEntry > &rasterEntries, const QgsCoordinateTransformContext &transformContext)
QgsRasterCalculator constructor.
QString lastError() const
Returns a description of the last error encountered.
Result processCalculation(QgsFeedback *feedback=nullptr)
Starts the calculation and writes a new raster.
Result
Result of the calculation.
@ InputLayerError
Error reading input layer.
@ ParserError
Error parsing formula.
@ CalculationError
Error occurred while performing calculation.
@ MemoryError
Error allocating memory for result.
@ Canceled
User canceled calculation.
@ Success
Calculation successful.
@ BandError
Invalid band number for input.
@ CreateOutputError
Error creating output data file.
QgsRasterBlock * block(int bandNo, const QgsRectangle &boundingBox, int width, int height, QgsRasterBlockFeedback *feedback=nullptr) override
Read block of data using given extent and size.
virtual Qgis::RasterInterfaceCapabilities capabilities() const
Returns the capabilities supported by the interface.
virtual bool setInput(QgsRasterInterface *input)
Set input.
Represents a raster layer.
int bandCount() const
Returns the number of bands in this layer.
QgsRasterDataProvider * dataProvider() override
Returns the source data provider.
Represents a matrix in a raster calculator operation.
bool isNumber() const
Returns true if matrix is 1x1 (=scalar number)
void setNodataValue(double d)
double * data()
Returns data array (but not ownership)
double number() const
Implements approximate projection support for optimised raster transformation.
void setPrecision(Precision precision)
QgsRasterBlock * block(int bandNo, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback=nullptr) override
Read block of data using given extent and size.
@ Exact
Exact, precise but slow.
Q_DECL_DEPRECATED void setCrs(const QgsCoordinateReferenceSystem &srcCRS, const QgsCoordinateReferenceSystem &destCRS, int srcDatumTransform=-1, int destDatumTransform=-1)
Sets the source and destination CRS.
A rectangle specified with double values.
double xMinimum() const
Returns the x minimum value (left side of rectangle).
void setYMinimum(double y)
Set the minimum y value.
double width() const
Returns the width of the rectangle.
double yMaximum() const
Returns the y maximum value (top side of rectangle).
void setYMaximum(double y)
Set the maximum y value.
double height() const
Returns the height of the rectangle.
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.
std::unique_ptr< std::remove_pointer< GDALDatasetH >::type, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
#define QgsDebugError(str)
Definition qgslogger.h:38
int CPL_STDCALL GdalProgressCallback(double dfComplete, const char *pszMessage, void *pProgressArg)
const QgsCoordinateReferenceSystem & outputCrs
const QString & typeName
Tiny smart-pointer-like wrapper around CPLMalloc and CPLFree: this is needed because OpenCL C++ API m...