QGIS API Documentation 3.28.0-Firenze (ed3ad0430f)
qgsalgorithmrasterstackposition.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsrasterstackposition.cpp
3 ---------------------
4 begin : July 2020
5 copyright : (C) 2020 by Clemens Raffler
6 email : clemens dot raffler 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 "qgsrasterprojector.h"
20#include "qgsrasterfilewriter.h"
22
24
25//
26//QgsRasterFrequencyByComparisonOperatorBase
27//
28
29QString QgsRasterStackPositionAlgorithmBase::group() const
30{
31 return QObject::tr( "Raster analysis" );
32}
33
34QString QgsRasterStackPositionAlgorithmBase::groupId() const
35{
36 return QStringLiteral( "rasteranalysis" );
37}
38
39void QgsRasterStackPositionAlgorithmBase::initAlgorithm( const QVariantMap & )
40{
41 addParameter( new QgsProcessingParameterMultipleLayers( QStringLiteral( "INPUT_RASTERS" ),
42 QObject::tr( "Input raster layers" ), QgsProcessing::TypeRaster ) );
43
44 addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "REFERENCE_LAYER" ), QObject::tr( "Reference layer" ) ) );
45
46 addParameter( new QgsProcessingParameterBoolean( QStringLiteral( "IGNORE_NODATA" ), QObject::tr( "Ignore NoData values" ), false ) );
47
48 std::unique_ptr< QgsProcessingParameterNumber > output_nodata_parameter = std::make_unique< QgsProcessingParameterNumber >( QStringLiteral( "OUTPUT_NODATA_VALUE" ), QObject::tr( "Output NoData value" ), QgsProcessingParameterNumber::Double, -9999, true );
49 output_nodata_parameter->setFlags( output_nodata_parameter->flags() | QgsProcessingParameterDefinition::FlagAdvanced );
50 addParameter( output_nodata_parameter.release() );
51
52 addParameter( new QgsProcessingParameterRasterDestination( QStringLiteral( "OUTPUT" ),
53 QObject::tr( "Output layer" ) ) );
54 addOutput( new QgsProcessingOutputString( QStringLiteral( "EXTENT" ), QObject::tr( "Extent" ) ) );
55 addOutput( new QgsProcessingOutputString( QStringLiteral( "CRS_AUTHID" ), QObject::tr( "CRS authority identifier" ) ) );
56 addOutput( new QgsProcessingOutputNumber( QStringLiteral( "WIDTH_IN_PIXELS" ), QObject::tr( "Width in pixels" ) ) );
57 addOutput( new QgsProcessingOutputNumber( QStringLiteral( "HEIGHT_IN_PIXELS" ), QObject::tr( "Height in pixels" ) ) );
58 addOutput( new QgsProcessingOutputNumber( QStringLiteral( "TOTAL_PIXEL_COUNT" ), QObject::tr( "Total pixel count" ) ) );
59}
60
61bool QgsRasterStackPositionAlgorithmBase::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
62{
63 QgsRasterLayer *referenceLayer = parameterAsRasterLayer( parameters, QStringLiteral( "REFERENCE_LAYER" ), context );
64 if ( !referenceLayer )
65 throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "REFERENCE_LAYER" ) ) );
66
67 mIgnoreNoData = parameterAsBool( parameters, QStringLiteral( "IGNORE_NODATA" ), context );
68 mNoDataValue = parameterAsDouble( parameters, QStringLiteral( "OUTPUT_NODATA_VALUE" ), context );
69
70 mCrs = referenceLayer->crs();
71 mRasterUnitsPerPixelX = referenceLayer->rasterUnitsPerPixelX();
72 mRasterUnitsPerPixelY = referenceLayer->rasterUnitsPerPixelY();
73 mLayerWidth = referenceLayer->width();
74 mLayerHeight = referenceLayer->height();
75 mExtent = referenceLayer->extent();
76
77 const QList< QgsMapLayer * > layers = parameterAsLayerList( parameters, QStringLiteral( "INPUT_RASTERS" ), context );
78 QList< QgsRasterLayer * > rasterLayers;
79 rasterLayers.reserve( layers.count() );
80 for ( QgsMapLayer *l : layers )
81 {
82 if ( feedback->isCanceled() )
83 break; //in case some slow data sources are loaded
84
85 if ( l->type() == QgsMapLayerType::RasterLayer )
86 {
87 QgsRasterLayer *layer = qobject_cast< QgsRasterLayer * >( l );
88 QgsRasterAnalysisUtils::RasterLogicInput input;
89 const int band = 1; //could be made dynamic
90 input.hasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( band );
91 input.sourceDataProvider.reset( layer->dataProvider()->clone() );
92 input.interface = input.sourceDataProvider.get();
93 // add projector if necessary
94 if ( layer->crs() != mCrs )
95 {
96 input.projector = std::make_unique< QgsRasterProjector >();
97 input.projector->setInput( input.sourceDataProvider.get() );
98 input.projector->setCrs( layer->crs(), mCrs, context.transformContext() );
99 input.interface = input.projector.get();
100 }
101 mInputs.emplace_back( std::move( input ) );
102 }
103 }
104
105 return true;
106}
107
108QVariantMap QgsRasterStackPositionAlgorithmBase::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
109{
110 const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context );
111 const QFileInfo fi( outputFile );
112 const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
113
114 std::unique_ptr< QgsRasterFileWriter > writer = std::make_unique< QgsRasterFileWriter >( outputFile );
115 writer->setOutputProviderKey( QStringLiteral( "gdal" ) );
116 writer->setOutputFormat( outputFormat );
117 std::unique_ptr<QgsRasterDataProvider > provider( writer->createOneBandRaster( Qgis::DataType::Int32, mLayerWidth, mLayerHeight, mExtent, mCrs ) );
118 if ( !provider )
119 throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( outputFile ) );
120 if ( !provider->isValid() )
121 throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( outputFile, provider->error().message( QgsErrorMessage::Text ) ) );
122
123 provider->setNoDataValue( 1, mNoDataValue );
124 const qgssize layerSize = static_cast< qgssize >( mLayerWidth ) * static_cast< qgssize >( mLayerHeight );
125
128 const int nbBlocksWidth = static_cast< int>( std::ceil( 1.0 * mLayerWidth / maxWidth ) );
129 const int nbBlocksHeight = static_cast< int >( std::ceil( 1.0 * mLayerHeight / maxHeight ) );
130 const int nbBlocks = nbBlocksWidth * nbBlocksHeight;
131 provider->setEditable( true );
132
133 QgsRasterIterator iter( provider.get() );
134 iter.startRasterRead( 1, mLayerWidth, mLayerHeight, mExtent );
135 int iterLeft = 0;
136 int iterTop = 0;
137 int iterCols = 0;
138 int iterRows = 0;
139 QgsRectangle blockExtent;
140
141 std::unique_ptr< QgsRasterBlock > outputBlock;
142 while ( iter.readNextRasterPart( 1, iterCols, iterRows, outputBlock, iterLeft, iterTop, &blockExtent ) )
143 {
144 std::vector< std::unique_ptr< QgsRasterBlock > > inputBlocks;
145 for ( const QgsRasterAnalysisUtils::RasterLogicInput &i : mInputs )
146 {
147 if ( feedback->isCanceled() )
148 break; //in case some slow data sources are loaded
149 for ( const int band : i.bands )
150 {
151 if ( feedback->isCanceled() )
152 break; //in case some slow data sources are loaded
153 std::unique_ptr< QgsRasterBlock > b( i.interface->block( band, blockExtent, iterCols, iterRows ) );
154 inputBlocks.emplace_back( std::move( b ) );
155 }
156 }
157
158 feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
159 for ( int row = 0; row < iterRows; row++ )
160 {
161 if ( feedback->isCanceled() )
162 break;
163
164 for ( int col = 0; col < iterCols; col++ )
165 {
166 bool noDataInStack = false;
167
168 if ( !inputBlocks.empty() )
169 {
170 const int position = findPosition( inputBlocks, row, col, noDataInStack );
171
172 if ( position == -1 || ( noDataInStack && !mIgnoreNoData ) )
173 {
174 //output cell will always be NoData if NoData occurs the current raster cell
175 //of the input blocks and NoData is not ignored
176 //this saves unnecessary iterations on the cellValueStack
177 outputBlock->setValue( row, col, mNoDataValue );
178 }
179 else
180 {
181 outputBlock->setValue( row, col, position );
182 }
183 }
184 else
185 {
186 outputBlock->setValue( row, col, mNoDataValue );
187 }
188 }
189 }
190 provider->writeBlock( outputBlock.get(), 1, iterLeft, iterTop );
191 }
192 provider->setEditable( false );
193
194 QVariantMap outputs;
195 outputs.insert( QStringLiteral( "EXTENT" ), mExtent.toString() );
196 outputs.insert( QStringLiteral( "CRS_AUTHID" ), mCrs.authid() );
197 outputs.insert( QStringLiteral( "WIDTH_IN_PIXELS" ), mLayerWidth );
198 outputs.insert( QStringLiteral( "HEIGHT_IN_PIXELS" ), mLayerHeight );
199 outputs.insert( QStringLiteral( "TOTAL_PIXEL_COUNT" ), layerSize );
200 outputs.insert( QStringLiteral( "OUTPUT" ), outputFile );
201
202 return outputs;
203}
204
205//
206// QgsRasterStackLowestPositionAlgorithm
207//
208QString QgsRasterStackLowestPositionAlgorithm::displayName() const
209{
210 return QObject::tr( "Lowest position in raster stack" );
211}
212
213QString QgsRasterStackLowestPositionAlgorithm::name() const
214{
215 return QStringLiteral( "lowestpositioninrasterstack" );
216}
217
218QStringList QgsRasterStackLowestPositionAlgorithm::tags() const
219{
220 return QObject::tr( "cell,lowest,position,pixel,stack" ).split( ',' );
221}
222
223QString QgsRasterStackLowestPositionAlgorithm::shortHelpString() const
224{
225 return QObject::tr( "The lowest position algorithm evaluates on a cell-by-cell basis the position "
226 "of the raster with the lowest value in a stack of rasters. Position counts start "
227 "with 1 and range to the total number of input rasters. The order of the input "
228 "rasters is relevant for the algorithm. If multiple rasters feature the lowest value, "
229 "the first raster will be used for the position value.\n "
230 "If multiband rasters are used in the data raster stack, the algorithm will always "
231 "perform the analysis on the first band of the rasters - use GDAL to use other bands in the analysis. "
232 "Any NoData cells in the raster layer stack will result in a NoData cell "
233 "in the output raster unless the \"ignore NoData\" parameter is checked. "
234 "The output NoData value can be set manually. The output rasters extent and resolution "
235 "is defined by a reference raster layer and is always of int32 type." );
236}
237
238QgsRasterStackLowestPositionAlgorithm *QgsRasterStackLowestPositionAlgorithm::createInstance() const
239{
240 return new QgsRasterStackLowestPositionAlgorithm();
241}
242
243int QgsRasterStackLowestPositionAlgorithm::findPosition( std::vector< std::unique_ptr<QgsRasterBlock> > &inputBlocks, int &row, int &col, bool &noDataInRasterBlockStack )
244{
245 int lowestPosition = 0;
246
247 //auxiliary variables
248 const int inputBlocksCount = inputBlocks.size();
249 int currentPosition = 0;
250 int noDataCount = 0;
251 double firstValue = mNoDataValue;
252 bool firstValueIsNoData = true;
253
254 while ( firstValueIsNoData && ( currentPosition < inputBlocksCount ) )
255 {
256 //check if all blocks are nodata/invalid
257 std::unique_ptr<QgsRasterBlock> &firstBlock = inputBlocks.at( currentPosition );
258 firstValue = firstBlock->valueAndNoData( row, col, firstValueIsNoData );
259
260 if ( !firstBlock->isValid() || firstValueIsNoData )
261 {
262 noDataInRasterBlockStack = true;
263 noDataCount++;
264 }
265 else
266 {
267 lowestPosition = currentPosition;
268 }
269 currentPosition++;
270 }
271
272 if ( noDataCount == inputBlocksCount )
273 {
274 noDataInRasterBlockStack = true;
275 return -1; //all blocks are NoData
276 }
277 else
278 {
279 //scan for the lowest value
280 while ( currentPosition < inputBlocksCount )
281 {
282 std::unique_ptr< QgsRasterBlock > &currentBlock = inputBlocks.at( currentPosition );
283
284 bool currentValueIsNoData = false;
285 const double currentValue = currentBlock->valueAndNoData( row, col, currentValueIsNoData );
286
287 if ( !currentBlock->isValid() || currentValueIsNoData )
288 {
289 noDataInRasterBlockStack = true;
290 noDataCount++;
291 }
292 else
293 {
294 if ( currentValue < firstValue )
295 {
296 firstValue = currentValue;
297 lowestPosition = currentPosition;
298 }
299 }
300 currentPosition++;
301 }
302 }
303 //the ArcGIS implementation uses 1 for first position value instead of 0 as in standard c++
304 return ++lowestPosition; //therefore ++
305}
306
307//
308// QgsRasterStackHighestPositionAlgorithmAlgorithm
309//
310
311QString QgsRasterStackHighestPositionAlgorithm::displayName() const
312{
313 return QObject::tr( "Highest position in raster stack" );
314}
315
316QString QgsRasterStackHighestPositionAlgorithm::name() const
317{
318 return QStringLiteral( "highestpositioninrasterstack" );
319}
320
321QStringList QgsRasterStackHighestPositionAlgorithm::tags() const
322{
323 return QObject::tr( "cell,highest,position,pixel,stack" ).split( ',' );
324}
325
326QString QgsRasterStackHighestPositionAlgorithm::shortHelpString() const
327{
328 return QObject::tr( "The highest position algorithm evaluates on a cell-by-cell basis the position "
329 "of the raster with the highest value in a stack of rasters. Position counts start "
330 "with 1 and range to the total number of input rasters. The order of the input "
331 "rasters is relevant for the algorithm. If multiple rasters feature the highest value, "
332 "the first raster will be used for the position value.\n "
333 "If multiband rasters are used in the data raster stack, the algorithm will always "
334 "perform the analysis on the first band of the rasters - use GDAL to use other bands in the analysis. "
335 "Any NoData cells in the raster layer stack will result in a NoData cell "
336 "in the output raster unless the \"ignore NoData\" parameter is checked. "
337 "The output NoData value can be set manually. The output rasters extent and resolution "
338 "is defined by a reference raster layer and is always of int32 type." );
339}
340
341QgsRasterStackHighestPositionAlgorithm *QgsRasterStackHighestPositionAlgorithm::createInstance() const
342{
343 return new QgsRasterStackHighestPositionAlgorithm();
344}
345
346int QgsRasterStackHighestPositionAlgorithm::findPosition( std::vector< std::unique_ptr< QgsRasterBlock> > &inputBlocks, int &row, int &col, bool &noDataInRasterBlockStack )
347{
348 int highestPosition = 0;
349
350 //auxiliary variables
351 const int inputBlocksCount = inputBlocks.size();
352 int currentPosition = 0;
353 int noDataCount = 0;
354 double firstValue = mNoDataValue;
355 bool firstValueIsNoData = true;
356
357 while ( firstValueIsNoData && ( currentPosition < inputBlocksCount ) )
358 {
359 //check if all blocks are nodata/invalid
360 std::unique_ptr<QgsRasterBlock> &firstBlock = inputBlocks.at( currentPosition );
361 firstValue = firstBlock->valueAndNoData( row, col, firstValueIsNoData );
362
363 if ( !firstBlock->isValid() || firstValueIsNoData )
364 {
365 noDataInRasterBlockStack = true;
366 noDataCount++;
367 }
368 else
369 {
370 highestPosition = currentPosition;
371 }
372
373 currentPosition++;
374 }
375
376 if ( noDataCount == inputBlocksCount )
377 {
378 noDataInRasterBlockStack = true;
379 return -1; //all blocks are NoData
380 }
381 else
382 {
383 //scan for the lowest value
384 while ( currentPosition < inputBlocksCount )
385 {
386 std::unique_ptr< QgsRasterBlock > &currentBlock = inputBlocks.at( currentPosition );
387
388 bool currentValueIsNoData = false;
389 const double currentValue = currentBlock->valueAndNoData( row, col, currentValueIsNoData );
390
391 if ( !currentBlock->isValid() || currentValueIsNoData )
392 {
393 noDataInRasterBlockStack = true;
394 noDataCount++;
395 }
396 else
397 {
398 if ( currentValue > firstValue )
399 {
400 firstValue = currentValue;
401 highestPosition = currentPosition;
402 }
403 }
404 currentPosition++;
405 }
406 }
407 //the ArcGIS implementation uses 1 for first position value instead of 0 as in standard c++
408 return ++highestPosition; //therefore ++
409}
410
412
@ Int32
Thirty two bit signed integer (qint32)
bool isCanceled() const SIP_HOLDGIL
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
Base class for all map layer types.
Definition: qgsmaplayer.h:73
virtual QgsRectangle extent() const
Returns the extent of the layer.
QgsCoordinateReferenceSystem crs
Definition: qgsmaplayer.h:79
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.
Definition: qgsexception.h:83
Base class for providing feedback from a processing algorithm.
A numeric output for processing algorithms.
A string output for processing algorithms.
A boolean parameter for processing algorithms.
@ FlagAdvanced
Parameter is an advanced parameter which should be hidden from users by default.
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 ...
A raster layer parameter for processing algorithms.
@ TypeRaster
Raster layers.
Definition: qgsprocessing.h:52
QgsRasterDataProvider * clone() const override=0
Clone itself, create deep copy.
virtual bool sourceHasNoDataValue(int bandNo) const
Returns true if source band has no data value.
static QString driverForExtension(const QString &extension)
Returns the GDAL driver name for a specified file extension.
Iterator for sequentially processing raster cells.
static const int DEFAULT_MAXIMUM_TILE_WIDTH
Default maximum tile width.
static const int DEFAULT_MAXIMUM_TILE_HEIGHT
Default maximum tile height.
Represents a raster layer.
int height() const
Returns the height of the (unclipped) raster.
double rasterUnitsPerPixelX() const
Returns the number of raster units per each raster pixel in X axis.
QgsRasterDataProvider * dataProvider() override
Returns the source data provider.
double rasterUnitsPerPixelY() const
Returns the number of raster units per each raster pixel in Y axis.
int width() const
Returns the width of the (unclipped) raster.
A rectangle specified with double values.
Definition: qgsrectangle.h:42
@ RasterLayer
Raster layer.
unsigned long long qgssize
Qgssize is used instead of size_t, because size_t is stdlib type, unknown by SIP, and it would be har...
Definition: qgis.h:3032