QGIS API Documentation 3.43.0-Master (e01d6d7c4c0)
qgsalgorithmrastercalculator.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsalgorithmrastercalculator.cpp
3 ---------------------
4 begin : July 2023
5 copyright : (C) 2023 by Alexander Bruy
6 email : alexander dot bruy at gmail dot com
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
19#include "qgsrasterfilewriter.h"
20#include "qgsrastercalculator.h"
21
23
24Qgis::ProcessingAlgorithmFlags QgsRasterCalculatorAlgorithm::flags() const
25{
27}
28
29QString QgsRasterCalculatorAlgorithm::name() const
30{
31 return QStringLiteral( "rastercalc" );
32}
33
34QString QgsRasterCalculatorAlgorithm::displayName() const
35{
36 return QObject::tr( "Raster calculator" );
37}
38
39QStringList QgsRasterCalculatorAlgorithm::tags() const
40{
41 return QObject::tr( "raster,calculator" ).split( ',' );
42}
43
44QString QgsRasterCalculatorAlgorithm::group() const
45{
46 return QObject::tr( "Raster analysis" );
47}
48
49QString QgsRasterCalculatorAlgorithm::groupId() const
50{
51 return QStringLiteral( "rasteranalysis" );
52}
53
54QString QgsRasterCalculatorAlgorithm::shortHelpString() const
55{
56 return QObject::tr( "This algorithm performs algebraic operations using raster layers." );
57}
58
59QString QgsRasterCalculatorAlgorithm::shortDescription() const
60{
61 return QObject::tr( "Performs algebraic operations using raster layers." );
62}
63
64QgsRasterCalculatorAlgorithm *QgsRasterCalculatorAlgorithm::createInstance() const
65{
66 return new QgsRasterCalculatorAlgorithm();
67}
68
69void QgsRasterCalculatorAlgorithm::initAlgorithm( const QVariantMap & )
70{
71 addParameter( new QgsProcessingParameterMultipleLayers( QStringLiteral( "LAYERS" ), QObject::tr( "Input layers" ), Qgis::ProcessingSourceType::Raster ) );
72 addParameter( new QgsProcessingParameterExpression( QStringLiteral( "EXPRESSION" ), QObject::tr( "Expression" ), QVariant(), QStringLiteral( "LAYERS" ), false, Qgis::ExpressionType::RasterCalculator ) );
73 auto extentParam = std::make_unique<QgsProcessingParameterExtent>( QStringLiteral( "EXTENT" ), QObject::tr( "Output extent" ), QVariant(), true );
74 extentParam->setHelp( QObject::tr( "Extent of the output layer. If not specified, the extent will be the overall extent of all input layers" ) );
75 addParameter( extentParam.release() );
76 auto cellSizeParam = std::make_unique<QgsProcessingParameterNumber>( QStringLiteral( "CELL_SIZE" ), QObject::tr( "Output cell size (leave empty to set automatically)" ), Qgis::ProcessingNumberParameterType::Double, QVariant(), true, 0.0 );
77 cellSizeParam->setHelp( QObject::tr( "Cell size of the output layer. If not specified, the smallest cell size from the input layers will be used" ) );
78 addParameter( cellSizeParam.release() );
79 auto crsParam = std::make_unique<QgsProcessingParameterCrs>( QStringLiteral( "CRS" ), QObject::tr( "Output CRS" ), QVariant(), true );
80 crsParam->setHelp( QObject::tr( "CRS of the output layer. If not specified, the CRS of the first input layer will be used" ) );
81 addParameter( crsParam.release() );
82
83 auto creationOptsParam = std::make_unique<QgsProcessingParameterString>( QStringLiteral( "CREATION_OPTIONS" ), QObject::tr( "Creation options" ), QVariant(), false, true );
84 creationOptsParam->setMetadata( QVariantMap( { { QStringLiteral( "widget_wrapper" ), QVariantMap( { { QStringLiteral( "widget_type" ), QStringLiteral( "rasteroptions" ) } } ) } } ) );
85 creationOptsParam->setFlags( creationOptsParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
86 addParameter( creationOptsParam.release() );
87
88 addParameter( new QgsProcessingParameterRasterDestination( QStringLiteral( "OUTPUT" ), QObject::tr( "Calculated" ) ) );
89}
90
91bool QgsRasterCalculatorAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
92{
93 const QList<QgsMapLayer *> layers = parameterAsLayerList( parameters, QStringLiteral( "LAYERS" ), context );
94
95 for ( const QgsMapLayer *layer : std::as_const( layers ) )
96 {
97 QgsMapLayer *clonedLayer { layer->clone() };
98 clonedLayer->moveToThread( nullptr );
99 mLayers << clonedLayer;
100 }
101
102 if ( mLayers.isEmpty() )
103 {
104 feedback->reportError( QObject::tr( "No layers selected" ), false );
105 return false;
106 }
107
108 return true;
109}
110
111QVariantMap QgsRasterCalculatorAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
112{
113 for ( QgsMapLayer *layer : std::as_const( mLayers ) )
114 {
115 layer->moveToThread( QThread::currentThread() );
116 }
117
119 if ( parameters.value( QStringLiteral( "CRS" ) ).isValid() )
120 {
121 crs = parameterAsCrs( parameters, QStringLiteral( "CRS" ), context );
122 }
123 else
124 {
125 crs = mLayers.at( 0 )->crs();
126 }
127
128 QgsRectangle bbox;
129 if ( parameters.value( QStringLiteral( "EXTENT" ) ).isValid() )
130 {
131 bbox = parameterAsExtent( parameters, QStringLiteral( "EXTENT" ), context, crs );
132 }
133 else
134 {
135 bbox = QgsProcessingUtils::combineLayerExtents( mLayers, crs, context );
136 }
137
138 double minCellSize = 1e9;
139
140 QVector<QgsRasterCalculatorEntry> entries;
141 for ( QgsMapLayer *layer : mLayers )
142 {
143 QgsRasterLayer *rLayer = qobject_cast<QgsRasterLayer *>( layer );
144 if ( !rLayer )
145 {
146 continue;
147 }
148
149 const int nBands = rLayer->dataProvider()->bandCount();
150 for ( int i = 0; i < nBands; ++i )
151 {
153 entry.ref = QStringLiteral( "%1@%2" ).arg( rLayer->name() ).arg( i + 1 );
154 entry.raster = rLayer;
155 entry.bandNumber = i + 1;
156 entries << entry;
157 }
158
159 QgsRectangle ext = rLayer->extent();
160 if ( rLayer->crs() != crs )
161 {
162 QgsCoordinateTransform ct( rLayer->crs(), crs, context.transformContext() );
163 ext = ct.transformBoundingBox( ext );
164 }
165
166 double cellSize = ( ext.xMaximum() - ext.xMinimum() ) / rLayer->width();
167 if ( cellSize < minCellSize )
168 {
169 minCellSize = cellSize;
170 }
171 }
172
173 double cellSize = parameterAsDouble( parameters, QStringLiteral( "CELL_SIZE" ), context );
174 if ( cellSize == 0 )
175 {
176 cellSize = minCellSize;
177 }
178
179 const QString creationOptions = parameterAsString( parameters, QStringLiteral( "CREATION_OPTIONS" ), context ).trimmed();
180 const QString expression = parameterAsExpression( parameters, QStringLiteral( "EXPRESSION" ), context );
181 const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context );
182 const QFileInfo fi( outputFile );
183 const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
184
185 double width = std::round( ( bbox.xMaximum() - bbox.xMinimum() ) / cellSize );
186 double height = std::round( ( bbox.yMaximum() - bbox.yMinimum() ) / cellSize );
187
188 QgsRasterCalculator calc( expression, outputFile, outputFormat, bbox, crs, width, height, entries, context.transformContext() );
189 calc.setCreationOptions( creationOptions.split( '|', Qt::SplitBehaviorFlags::SkipEmptyParts ) );
190 QgsRasterCalculator::Result result = calc.processCalculation( feedback );
191 qDeleteAll( mLayers );
192 mLayers.clear();
193 switch ( result )
194 {
196 throw QgsProcessingException( QObject::tr( "Error creating output file." ) );
198 throw QgsProcessingException( QObject::tr( "Error reading input layer." ) );
200 throw QgsProcessingException( QObject::tr( "Error parsing formula." ) );
202 throw QgsProcessingException( QObject::tr( "Error allocating memory for result." ) );
204 throw QgsProcessingException( QObject::tr( "Invalid band number for input." ) );
206 throw QgsProcessingException( QObject::tr( "Error occurred while performing calculation." ) );
207 default:
208 break;
209 }
210
211 QVariantMap outputs;
212 outputs.insert( QStringLiteral( "OUTPUT" ), outputFile );
213 return outputs;
214}
215
216Qgis::ProcessingAlgorithmFlags QgsRasterCalculatorModelerAlgorithm::flags() const
217{
219}
220
221QString QgsRasterCalculatorModelerAlgorithm::name() const
222{
223 return QStringLiteral( "modelerrastercalc" );
224}
225
226QString QgsRasterCalculatorModelerAlgorithm::displayName() const
227{
228 return QObject::tr( "Raster calculator" );
229}
230
231QStringList QgsRasterCalculatorModelerAlgorithm::tags() const
232{
233 return QObject::tr( "raster,calculator" ).split( ',' );
234}
235
236QString QgsRasterCalculatorModelerAlgorithm::group() const
237{
238 return QObject::tr( "Raster analysis" );
239}
240
241QString QgsRasterCalculatorModelerAlgorithm::groupId() const
242{
243 return QStringLiteral( "rasteranalysis" );
244}
245
246QgsRasterCalculatorModelerAlgorithm *QgsRasterCalculatorModelerAlgorithm::createInstance() const
247{
248 return new QgsRasterCalculatorModelerAlgorithm();
249}
250
251QVariantMap QgsRasterCalculatorModelerAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
252{
253 for ( QgsMapLayer *layer : std::as_const( mLayers ) )
254 {
255 layer->moveToThread( QThread::currentThread() );
256 }
257
259 if ( parameters.value( QStringLiteral( "CRS" ) ).isValid() )
260 {
261 crs = parameterAsCrs( parameters, QStringLiteral( "CRS" ), context );
262 }
263 else
264 {
265 crs = mLayers.at( 0 )->crs();
266 }
267
268 QgsRectangle bbox;
269 if ( parameters.value( QStringLiteral( "EXTENT" ) ).isValid() )
270 {
271 bbox = parameterAsExtent( parameters, QStringLiteral( "EXTENT" ), context, crs );
272 }
273 else
274 {
275 bbox = QgsProcessingUtils::combineLayerExtents( mLayers, crs, context );
276 }
277
278 double minCellSize = 1e9;
279
280 QVector<QgsRasterCalculatorEntry> entries;
281 int n = 0;
282 for ( QgsMapLayer *layer : mLayers )
283 {
284 QgsRasterLayer *rLayer = qobject_cast<QgsRasterLayer *>( layer );
285 if ( !rLayer )
286 {
287 continue;
288 }
289
290 n++;
291 const int nBands = rLayer->dataProvider()->bandCount();
292 for ( int i = 0; i < nBands; ++i )
293 {
295 entry.ref = QStringLiteral( "%1@%2" ).arg( indexToName( n ) ).arg( i + 1 );
296 entry.raster = rLayer;
297 entry.bandNumber = i + 1;
298 entries << entry;
299 }
300
301 QgsRectangle ext = rLayer->extent();
302 if ( rLayer->crs() != crs )
303 {
304 QgsCoordinateTransform ct( rLayer->crs(), crs, context.transformContext() );
305 ext = ct.transformBoundingBox( ext );
306 }
307
308 double cellSize = ( ext.xMaximum() - ext.xMinimum() ) / rLayer->width();
309 if ( cellSize < minCellSize )
310 {
311 minCellSize = cellSize;
312 }
313 }
314
315 double cellSize = parameterAsDouble( parameters, QStringLiteral( "CELL_SIZE" ), context );
316 if ( cellSize == 0 )
317 {
318 cellSize = minCellSize;
319 }
320
321 const QString expression = parameterAsExpression( parameters, QStringLiteral( "EXPRESSION" ), context );
322 const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context );
323 const QFileInfo fi( outputFile );
324 const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
325
326 double width = std::round( ( bbox.xMaximum() - bbox.xMinimum() ) / cellSize );
327 double height = std::round( ( bbox.yMaximum() - bbox.yMinimum() ) / cellSize );
328
329 QgsRasterCalculator calc( expression, outputFile, outputFormat, bbox, crs, width, height, entries, context.transformContext() );
330 QgsRasterCalculator::Result result = calc.processCalculation( feedback );
331 qDeleteAll( mLayers );
332 mLayers.clear();
333 switch ( result )
334 {
336 throw QgsProcessingException( QObject::tr( "Error creating output file." ) );
338 throw QgsProcessingException( QObject::tr( "Error reading input layer." ) );
340 throw QgsProcessingException( QObject::tr( "Error parsing formula." ) );
342 throw QgsProcessingException( QObject::tr( "Error allocating memory for result." ) );
344 throw QgsProcessingException( QObject::tr( "Invalid band number for input." ) );
346 throw QgsProcessingException( QObject::tr( "Error occurred while performing calculation." ) );
347 default:
348 break;
349 }
350
351 QVariantMap outputs;
352 outputs.insert( QStringLiteral( "OUTPUT" ), outputFile );
353 return outputs;
354}
355
356QString QgsRasterCalculatorModelerAlgorithm::indexToName( int index ) const
357{
358 QString name;
359 int div = index;
360 int mod = 0;
361
362 while ( div > 0 )
363 {
364 mod = ( div - 1 ) % 26;
365 name = static_cast<char>( 65 + mod ) + name;
366 div = ( int ) ( ( div - mod ) / 26 );
367 }
368 return name;
369}
370
@ RasterCalculator
Raster calculator expression.
QFlags< ProcessingAlgorithmFlag > ProcessingAlgorithmFlags
Flags indicating how and when an algorithm operates and should be exposed to users.
Definition qgis.h:3476
@ HideFromToolbox
Algorithm should be hidden from the toolbox.
@ HideFromModeler
Algorithm should be hidden from the modeler.
@ Advanced
Parameter is an advanced parameter which should be hidden from users by default.
Represents a coordinate reference system (CRS).
Handles coordinate transforms between two coordinate systems.
Base class for all map layer types.
Definition qgsmaplayer.h:77
QString name
Definition qgsmaplayer.h:81
virtual QgsRectangle extent() const
Returns the extent of the layer.
QgsCoordinateReferenceSystem crs
Definition qgsmaplayer.h:84
virtual QgsMapLayer * clone() const =0
Returns a new instance equivalent to this one except for the id which is still unique.
virtual Qgis::ProcessingAlgorithmFlags flags() const
Returns the flags indicating how and when the algorithm operates and should be exposed to users.
Contains information about the context in which a processing algorithm is executed.
QgsCoordinateTransformContext transformContext() const
Returns the coordinate transform context.
Custom exception class for processing related exceptions.
Base class for providing feedback from a processing algorithm.
virtual void reportError(const QString &error, bool fatalError=false)
Reports that the algorithm encountered an error while executing.
An expression parameter for processing algorithms.
A parameter for processing algorithms which accepts multiple map layers.
A raster layer destination parameter, for specifying the destination path for a raster layer created ...
static QgsRectangle combineLayerExtents(const QList< QgsMapLayer * > &layers, const QgsCoordinateReferenceSystem &crs, QgsProcessingContext &context)
Combines the extent of several map layers.
Represents an individual raster layer/band number entry within a raster calculation.
QgsRasterLayer * raster
Raster layer associated with entry.
int bandNumber
Band number for entry.
QString ref
Name of entry.
Performs raster layer calculations.
Result
Result of the calculation.
@ InputLayerError
Error reading input layer.
@ BandError
Invalid band number for input.
@ CreateOutputError
Error creating output data file.
@ ParserError
Error parsing formula.
@ CalculationError
Error occurred while performing calculation.
@ MemoryError
Error allocating memory for result.
static QString driverForExtension(const QString &extension)
Returns the GDAL driver name for a specified file extension.
virtual int bandCount() const =0
Gets number of bands.
Represents a raster layer.
QgsRasterDataProvider * dataProvider() override
Returns the source data provider.
int width() const
Returns the width of the (unclipped) raster.
A rectangle specified with double values.
double xMinimum
double yMinimum
double xMaximum
double yMaximum
const QgsCoordinateReferenceSystem & crs