QGIS API Documentation  3.26.3-Buenos Aires (65e4edfdad)
qgsalgorithmrasterfrequencybycomparisonoperator.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsalgorithmrasterfrequencybycomparisonoperator.cpp
3  ---------------------
4  begin : June 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 QgsRasterFrequencyByComparisonOperatorBase::group() const
30 {
31  return QObject::tr( "Raster analysis" );
32 }
33 
34 QString QgsRasterFrequencyByComparisonOperatorBase::groupId() const
35 {
36  return QStringLiteral( "rasteranalysis" );
37 }
38 
39 void QgsRasterFrequencyByComparisonOperatorBase::initAlgorithm( const QVariantMap & )
40 {
41  addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT_VALUE_RASTER" ), QObject::tr( "Input value raster" ) ) );
42  addParameter( new QgsProcessingParameterBand( QStringLiteral( "INPUT_VALUE_RASTER_BAND" ), QObject::tr( "Value raster band" ), 1, QStringLiteral( "INPUT_VALUE_RASTER" ) ) );
43 
44 
45  addParameter( new QgsProcessingParameterMultipleLayers( QStringLiteral( "INPUT_RASTERS" ),
46  QObject::tr( "Input raster layers" ), QgsProcessing::TypeRaster ) );
47 
48  addParameter( new QgsProcessingParameterBoolean( QStringLiteral( "IGNORE_NODATA" ), QObject::tr( "Ignore NoData values" ), false ) );
49 
50  std::unique_ptr< QgsProcessingParameterNumber > output_nodata_parameter = std::make_unique< QgsProcessingParameterNumber >( QStringLiteral( "OUTPUT_NODATA_VALUE" ), QObject::tr( "Output NoData value" ), QgsProcessingParameterNumber::Double, -9999, true );
51  output_nodata_parameter->setFlags( output_nodata_parameter->flags() | QgsProcessingParameterDefinition::FlagAdvanced );
52  addParameter( output_nodata_parameter.release() );
53 
54  addParameter( new QgsProcessingParameterRasterDestination( QStringLiteral( "OUTPUT" ),
55  QObject::tr( "Output layer" ) ) );
56  addOutput( new QgsProcessingOutputNumber( QStringLiteral( "OCCURRENCE_COUNT" ), QObject::tr( "Count of value occurrences" ) ) );
57  addOutput( new QgsProcessingOutputNumber( QStringLiteral( "FOUND_LOCATIONS_COUNT" ), QObject::tr( "Count of cells with equal value occurrences" ) ) );
58  addOutput( new QgsProcessingOutputNumber( QStringLiteral( "MEAN_FREQUENCY_PER_LOCATION" ), QObject::tr( "Mean frequency at valid cell locations" ) ) );
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 
66 bool QgsRasterFrequencyByComparisonOperatorBase::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
67 {
68  QgsRasterLayer *inputValueRaster = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT_VALUE_RASTER" ), context );
69  if ( !inputValueRaster )
70  throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT_VALUE_RASTER" ) ) );
71 
72  mInputValueRasterBand = parameterAsInt( parameters, QStringLiteral( "INPUT_VALUE_RASTER_BAND" ), context );
73  mIgnoreNoData = parameterAsBool( parameters, QStringLiteral( "IGNORE_NODATA" ), context );
74 
75  mInputValueRasterInterface.reset( inputValueRaster->dataProvider()->clone() );
76  mNoDataValue = parameterAsDouble( parameters, QStringLiteral( "OUTPUT_NODATA_VALUE" ), context );
77  mCrs = inputValueRaster->crs();
78  mRasterUnitsPerPixelX = inputValueRaster->rasterUnitsPerPixelX();
79  mRasterUnitsPerPixelY = inputValueRaster->rasterUnitsPerPixelY();
80  mLayerWidth = inputValueRaster->width();
81  mLayerHeight = inputValueRaster->height();
82  mExtent = inputValueRaster->extent();
83 
84  const QList< QgsMapLayer * > layers = parameterAsLayerList( parameters, QStringLiteral( "INPUT_RASTERS" ), context );
85  QList< QgsRasterLayer * > rasterLayers;
86  rasterLayers.reserve( layers.count() );
87  for ( QgsMapLayer *l : layers )
88  {
89  if ( feedback->isCanceled() )
90  break; //in case some slow data sources are loaded
91 
92  if ( l->type() == QgsMapLayerType::RasterLayer )
93  {
94  QgsRasterLayer *layer = qobject_cast< QgsRasterLayer * >( l );
95  QgsRasterAnalysisUtils::RasterLogicInput input;
96  const int band = 1; //could be made dynamic
97  input.hasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( band );
98  input.sourceDataProvider.reset( layer->dataProvider()->clone() );
99  input.interface = input.sourceDataProvider.get();
100  // add projector if necessary
101  if ( layer->crs() != mCrs )
102  {
103  input.projector = std::make_unique< QgsRasterProjector >();
104  input.projector->setInput( input.sourceDataProvider.get() );
105  input.projector->setCrs( layer->crs(), mCrs, context.transformContext() );
106  input.interface = input.projector.get();
107  }
108  mInputs.emplace_back( std::move( input ) );
109  }
110  }
111 
112  return true;
113 }
114 
115 QVariantMap QgsRasterFrequencyByComparisonOperatorBase::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
116 {
117  const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context );
118  const QFileInfo fi( outputFile );
119  const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
120 
121  std::unique_ptr< QgsRasterFileWriter > writer = std::make_unique< QgsRasterFileWriter >( outputFile );
122  writer->setOutputProviderKey( QStringLiteral( "gdal" ) );
123  writer->setOutputFormat( outputFormat );
124  std::unique_ptr<QgsRasterDataProvider > provider( writer->createOneBandRaster( Qgis::DataType::Int32, mLayerWidth, mLayerHeight, mExtent, mCrs ) );
125  if ( !provider )
126  throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( outputFile ) );
127  if ( !provider->isValid() )
128  throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( outputFile, provider->error().message( QgsErrorMessage::Text ) ) );
129 
130  provider->setNoDataValue( 1, mNoDataValue );
131  const qgssize layerSize = static_cast< qgssize >( mLayerWidth ) * static_cast< qgssize >( mLayerHeight );
132 
135  const int nbBlocksWidth = static_cast< int>( std::ceil( 1.0 * mLayerWidth / maxWidth ) );
136  const int nbBlocksHeight = static_cast< int >( std::ceil( 1.0 * mLayerHeight / maxHeight ) );
137  const int nbBlocks = nbBlocksWidth * nbBlocksHeight;
138  provider->setEditable( true );
139 
140  QgsRasterIterator iter( mInputValueRasterInterface.get() );
141  iter.startRasterRead( mInputValueRasterBand, mLayerWidth, mLayerHeight, mExtent );
142  int iterLeft = 0;
143  int iterTop = 0;
144  int iterCols = 0;
145  int iterRows = 0;
146  QgsRectangle blockExtent;
147 
148  unsigned long long occurrenceCount = 0;
149  unsigned long long noDataLocationsCount = 0;
150  std::unique_ptr< QgsRasterBlock > inputBlock;
151  while ( iter.readNextRasterPart( 1, iterCols, iterRows, inputBlock, iterLeft, iterTop, &blockExtent ) )
152  {
153  std::vector< std::unique_ptr< QgsRasterBlock > > inputBlocks;
154  for ( const QgsRasterAnalysisUtils::RasterLogicInput &i : mInputs )
155  {
156  if ( feedback->isCanceled() )
157  break; //in case some slow data sources are loaded
158  for ( const int band : i.bands )
159  {
160  if ( feedback->isCanceled() )
161  break; //in case some slow data sources are loaded
162  std::unique_ptr< QgsRasterBlock > b( i.interface->block( band, blockExtent, iterCols, iterRows ) );
163  inputBlocks.emplace_back( std::move( b ) );
164  }
165  }
166 
167  std::unique_ptr< QgsRasterBlock > outputBlock = std::make_unique<QgsRasterBlock>( Qgis::DataType::Int32, iterCols, iterRows );
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 valueRasterCellIsNoData = false;
177  const double value = inputBlock->valueAndNoData( row, col, valueRasterCellIsNoData );
178 
179  if ( valueRasterCellIsNoData && !mIgnoreNoData )
180  {
181  //output cell will always be NoData if NoData occurs in valueRaster or cellValueStack and NoData is not ignored
182  //this saves unnecessary iterations on the cellValueStack
183  outputBlock->setValue( row, col, mNoDataValue );
184  noDataLocationsCount++;
185  }
186  else
187  {
188  bool noDataInStack = false;
189  const std::vector<double> cellValues = QgsRasterAnalysisUtils::getCellValuesFromBlockStack( inputBlocks, row, col, noDataInStack );
190 
191  if ( noDataInStack && !mIgnoreNoData )
192  {
193  outputBlock->setValue( row, col, mNoDataValue );
194  noDataLocationsCount++;
195  }
196  else
197  {
198  const int frequency = applyComparisonOperator( value, cellValues );
199  outputBlock->setValue( row, col, frequency );
200  occurrenceCount += frequency;
201  }
202  }
203  }
204  }
205  provider->writeBlock( outputBlock.get(), 1, iterLeft, iterTop );
206  }
207  provider->setEditable( false );
208 
209  const unsigned long long foundLocationsCount = layerSize - noDataLocationsCount;
210  const double meanEqualCountPerValidLocation = static_cast<double>( occurrenceCount ) / static_cast<double>( foundLocationsCount * mInputs.size() );
211 
212  QVariantMap outputs;
213  outputs.insert( QStringLiteral( "OCCURRENCE_COUNT" ), occurrenceCount );
214  outputs.insert( QStringLiteral( "FOUND_LOCATIONS_COUNT" ), foundLocationsCount );
215  outputs.insert( QStringLiteral( "MEAN_FREQUENCY_PER_LOCATION" ), meanEqualCountPerValidLocation );
216  outputs.insert( QStringLiteral( "EXTENT" ), mExtent.toString() );
217  outputs.insert( QStringLiteral( "CRS_AUTHID" ), mCrs.authid() );
218  outputs.insert( QStringLiteral( "WIDTH_IN_PIXELS" ), mLayerWidth );
219  outputs.insert( QStringLiteral( "HEIGHT_IN_PIXELS" ), mLayerHeight );
220  outputs.insert( QStringLiteral( "TOTAL_PIXEL_COUNT" ), layerSize );
221  outputs.insert( QStringLiteral( "OUTPUT" ), outputFile );
222 
223  return outputs;
224 }
225 
226 //
227 // QgsRasterFrequencyByEqualOperatorAlgorithm
228 //
229 
230 QString QgsRasterFrequencyByEqualOperatorAlgorithm::displayName() const
231 {
232  return QObject::tr( "Equal to frequency" );
233 }
234 
235 QString QgsRasterFrequencyByEqualOperatorAlgorithm::name() const
236 {
237  return QStringLiteral( "equaltofrequency" );
238 }
239 
240 QStringList QgsRasterFrequencyByEqualOperatorAlgorithm::tags() const
241 {
242  return QObject::tr( "cell,equal,frequency,pixel,stack" ).split( ',' );
243 }
244 
245 QString QgsRasterFrequencyByEqualOperatorAlgorithm::shortHelpString() const
246 {
247  return QObject::tr( "The Equal to frequency algorithm evaluates on a cell-by-cell basis the frequency "
248  "(number of times) the values of an input stack of rasters are equal "
249  "to the value of a value raster. \n "
250  "If multiband rasters are used in the data raster stack, the algorithm will always "
251  "perform the analysis on the first band of the rasters - use GDAL to use other bands in the analysis. "
252  "The input value layer serves as reference layer for the sample layers. "
253  "Any NoData cells in the value raster or the data layer stack will result in a NoData cell "
254  "in the output raster if the ignore NoData parameter is not checked. "
255  "The output NoData value can be set manually. The output rasters extent and resolution "
256  "is defined by the input raster layer and is always of int32 type." );
257 }
258 
259 QgsRasterFrequencyByEqualOperatorAlgorithm *QgsRasterFrequencyByEqualOperatorAlgorithm::createInstance() const
260 {
261  return new QgsRasterFrequencyByEqualOperatorAlgorithm();
262 }
263 
264 int QgsRasterFrequencyByEqualOperatorAlgorithm::applyComparisonOperator( double searchValue, std::vector<double> cellValueStack )
265 {
266  return static_cast<int>( std::count( cellValueStack.begin(), cellValueStack.end(), searchValue ) );
267 }
268 
269 //
270 // QgsRasterFrequencyByGreaterThanOperatorAlgorithm
271 //
272 
273 QString QgsRasterFrequencyByGreaterThanOperatorAlgorithm::displayName() const
274 {
275  return QObject::tr( "Greater than frequency" );
276 }
277 
278 QString QgsRasterFrequencyByGreaterThanOperatorAlgorithm::name() const
279 {
280  return QStringLiteral( "greaterthanfrequency" );
281 }
282 
283 QStringList QgsRasterFrequencyByGreaterThanOperatorAlgorithm::tags() const
284 {
285  return QObject::tr( "cell,greater,frequency,pixel,stack" ).split( ',' );
286 }
287 
288 QString QgsRasterFrequencyByGreaterThanOperatorAlgorithm::shortHelpString() const
289 {
290  return QObject::tr( "The Greater than frequency algorithm evaluates on a cell-by-cell basis the frequency "
291  "(number of times) the values of an input stack of rasters are greater than "
292  "the value of a value raster. \n "
293  "If multiband rasters are used in the data raster stack, the algorithm will always "
294  "perform the analysis on the first band of the rasters - use GDAL to use other bands in the analysis. "
295  "The input value layer serves as reference layer for the sample layers. "
296  "Any NoData cells in the value raster or the data layer stack will result in a NoData cell "
297  "in the output raster if the ignore NoData parameter is not checked. "
298  "The output NoData value can be set manually. The output rasters extent and resolution "
299  "is defined by the input raster layer and is always of int32 type." );
300 }
301 
302 QgsRasterFrequencyByGreaterThanOperatorAlgorithm *QgsRasterFrequencyByGreaterThanOperatorAlgorithm::createInstance() const
303 {
304  return new QgsRasterFrequencyByGreaterThanOperatorAlgorithm();
305 }
306 
307 int QgsRasterFrequencyByGreaterThanOperatorAlgorithm::applyComparisonOperator( double searchValue, std::vector<double> cellValueStack )
308 {
309  return static_cast<int>( std::count_if( cellValueStack.begin(), cellValueStack.end(), [&]( double const & stackValue ) { return stackValue > searchValue; } ) );
310 }
311 
312 //
313 // QgsRasterFrequencyByLessThanOperatorAlgorithm
314 //
315 
316 QString QgsRasterFrequencyByLessThanOperatorAlgorithm::displayName() const
317 {
318  return QObject::tr( "Less than frequency" );
319 }
320 
321 QString QgsRasterFrequencyByLessThanOperatorAlgorithm::name() const
322 {
323  return QStringLiteral( "lessthanfrequency" );
324 }
325 
326 QStringList QgsRasterFrequencyByLessThanOperatorAlgorithm::tags() const
327 {
328  return QObject::tr( "cell,less,lower,frequency,pixel,stack" ).split( ',' );
329 }
330 
331 QString QgsRasterFrequencyByLessThanOperatorAlgorithm::shortHelpString() const
332 {
333  return QObject::tr( "The Less than frequency algorithm evaluates on a cell-by-cell basis the frequency "
334  "(number of times) the values of an input stack of rasters are less than "
335  "the value of a value raster. \n "
336  "If multiband rasters are used in the data raster stack, the algorithm will always "
337  "perform the analysis on the first band of the rasters - use GDAL to use other bands in the analysis. "
338  "The input value layer serves as reference layer for the sample layers. "
339  "Any NoData cells in the value raster or the data layer stack will result in a NoData cell "
340  "in the output raster if the ignore NoData parameter is not checked. "
341  "The output NoData value can be set manually. The output rasters extent and resolution "
342  "is defined by the input raster layer and is always of int32 type." );
343 }
344 
345 QgsRasterFrequencyByLessThanOperatorAlgorithm *QgsRasterFrequencyByLessThanOperatorAlgorithm::createInstance() const
346 {
347  return new QgsRasterFrequencyByLessThanOperatorAlgorithm();
348 }
349 
350 int QgsRasterFrequencyByLessThanOperatorAlgorithm::applyComparisonOperator( double searchValue, std::vector<double> cellValueStack )
351 {
352  return static_cast<int>( std::count_if( cellValueStack.begin(), cellValueStack.end(), [&]( double const & stackValue ) { return stackValue < searchValue; } ) );
353 }
354 
356 
357 
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
qgsalgorithmrasterfrequencybycomparisonoperator.h
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
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
QgsProcessingParameterBand
A raster band parameter for Processing algorithms.
Definition: qgsprocessingparameters.h:3548
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