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