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