QGIS API Documentation  3.26.3-Buenos Aires (65e4edfdad)
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"
21 #include "qgsrasteranalysisutils.h"
22 
24 
25 //
26 //QgsRasterFrequencyByComparisonOperatorBase
27 //
28 
29 QString QgsRasterStackPositionAlgorithmBase::group() const
30 {
31  return QObject::tr( "Raster analysis" );
32 }
33 
34 QString QgsRasterStackPositionAlgorithmBase::groupId() const
35 {
36  return QStringLiteral( "rasteranalysis" );
37 }
38 
39 void 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 
61 bool 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 
108 QVariantMap 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 //
208 QString QgsRasterStackLowestPositionAlgorithm::displayName() const
209 {
210  return QObject::tr( "Lowest position in raster stack" );
211 }
212 
213 QString QgsRasterStackLowestPositionAlgorithm::name() const
214 {
215  return QStringLiteral( "lowestpositioninrasterstack" );
216 }
217 
218 QStringList QgsRasterStackLowestPositionAlgorithm::tags() const
219 {
220  return QObject::tr( "cell,lowest,position,pixel,stack" ).split( ',' );
221 }
222 
223 QString 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 
238 QgsRasterStackLowestPositionAlgorithm *QgsRasterStackLowestPositionAlgorithm::createInstance() const
239 {
240  return new QgsRasterStackLowestPositionAlgorithm();
241 }
242 
243 int 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 
311 QString QgsRasterStackHighestPositionAlgorithm::displayName() const
312 {
313  return QObject::tr( "Highest position in raster stack" );
314 }
315 
316 QString QgsRasterStackHighestPositionAlgorithm::name() const
317 {
318  return QStringLiteral( "highestpositioninrasterstack" );
319 }
320 
321 QStringList QgsRasterStackHighestPositionAlgorithm::tags() const
322 {
323  return QObject::tr( "cell,highest,position,pixel,stack" ).split( ',' );
324 }
325 
326 QString 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 
341 QgsRasterStackHighestPositionAlgorithm *QgsRasterStackHighestPositionAlgorithm::createInstance() const
342 {
343  return new QgsRasterStackHighestPositionAlgorithm();
344 }
345 
346 int 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 
QgsMapLayer::crs
QgsCoordinateReferenceSystem crs
Definition: qgsmaplayer.h:79
QgsFeedback::setProgress
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:76
qgsrasterprojector.h
QgsProcessingParameterNumber::Double
@ Double
Double/float values.
Definition: qgsprocessingparameters.h:2187
QgsProcessingFeedback
Base class for providing feedback from a processing algorithm.
Definition: qgsprocessingfeedback.h:37
QgsRasterFileWriter::driverForExtension
static QString driverForExtension(const QString &extension)
Returns the GDAL driver name for a specified file extension.
Definition: qgsrasterfilewriter.cpp:1066
QgsProcessingParameterRasterDestination
A raster layer destination parameter, for specifying the destination path for a raster layer created ...
Definition: qgsprocessingparameters.h:3390
QgsFeedback::isCanceled
bool isCanceled() const SIP_HOLDGIL
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:67
QgsProcessingParameterDefinition::FlagAdvanced
@ FlagAdvanced
Parameter is an advanced parameter which should be hidden from users by default.
Definition: qgsprocessingparameters.h:451
Qgis::DataType::Int32
@ Int32
Thirty two bit signed integer (qint32)
qgsrasteranalysisutils.h
QgsProcessingParameterMultipleLayers
A parameter for processing algorithms which accepts multiple map layers.
Definition: qgsprocessingparameters.h:2097
QgsRectangle
A rectangle specified with double values.
Definition: qgsrectangle.h:41
QgsProcessingOutputNumber
A numeric output for processing algorithms.
Definition: qgsprocessingoutputs.h:312
qgsrasterfilewriter.h
QgsRasterDataProvider::clone
QgsRasterDataProvider * clone() const override=0
Clone itself, create deep copy.
QgsRasterDataProvider::sourceHasNoDataValue
virtual bool sourceHasNoDataValue(int bandNo) const
Returns true if source band has no data value.
Definition: qgsrasterdataprovider.h:241
QgsRasterLayer::width
int width() const
Returns the width of the (unclipped) raster.
Definition: qgsrasterlayer.cpp:2618
QgsProcessingContext
Contains information about the context in which a processing algorithm is executed.
Definition: qgsprocessingcontext.h:46
QgsRasterLayer::height
int height() const
Returns the height of the (unclipped) raster.
Definition: qgsrasterlayer.cpp:2624
QgsMapLayer::extent
virtual QgsRectangle extent() const
Returns the extent of the layer.
Definition: qgsmaplayer.cpp:305
qgsalgorithmrasterstackposition.h
QgsProcessing::TypeRaster
@ TypeRaster
Raster layers.
Definition: qgsprocessing.h:52
QgsMapLayerType::RasterLayer
@ RasterLayer
Raster layer.
QgsRasterIterator::DEFAULT_MAXIMUM_TILE_HEIGHT
static const int DEFAULT_MAXIMUM_TILE_HEIGHT
Default maximum tile height.
Definition: qgsrasteriterator.h:170
QgsProcessingParameterRasterLayer
A raster layer parameter for processing algorithms.
Definition: qgsprocessingparameters.h:2495
QgsProcessingOutputString
A string output for processing algorithms.
Definition: qgsprocessingoutputs.h:334
QgsProcessingContext::transformContext
QgsCoordinateTransformContext transformContext() const
Returns the coordinate transform context.
Definition: qgsprocessingcontext.h:165
QgsRasterLayer
Represents a raster layer.
Definition: qgsrasterlayer.h:76
QgsRasterIterator
Iterator for sequentially processing raster cells.
Definition: qgsrasteriterator.h:34
QgsRasterLayer::rasterUnitsPerPixelY
double rasterUnitsPerPixelY() const
Returns the number of raster units per each raster pixel in Y axis.
Definition: qgsrasterlayer.cpp:582
QgsProcessingParameterBoolean
A boolean parameter for processing algorithms.
Definition: qgsprocessingparameters.h:1709
QgsMapLayer
Base class for all map layer types. This is the base class for all map layer types (vector,...
Definition: qgsmaplayer.h:72
QgsRasterIterator::DEFAULT_MAXIMUM_TILE_WIDTH
static const int DEFAULT_MAXIMUM_TILE_WIDTH
Default maximum tile width.
Definition: qgsrasteriterator.h:167
QgsErrorMessage::Text
@ Text
Definition: qgserror.h:38
QgsProcessingException
Custom exception class for processing related exceptions.
Definition: qgsexception.h:82
QgsRasterLayer::dataProvider
QgsRasterDataProvider * dataProvider() override
Returns the source data provider.
Definition: qgsrasterlayer.cpp:257
QgsRasterLayer::rasterUnitsPerPixelX
double rasterUnitsPerPixelX() const
Returns the number of raster units per each raster pixel in X axis.
Definition: qgsrasterlayer.cpp:566
qgssize
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:2791