QGIS API Documentation 3.28.0-Firenze (ed3ad0430f)
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"
22
24
25//
26//QgsRasterFrequencyByComparisonOperatorBase
27//
28
29QString QgsRasterFrequencyByComparisonOperatorBase::group() const
30{
31 return QObject::tr( "Raster analysis" );
32}
33
34QString QgsRasterFrequencyByComparisonOperatorBase::groupId() const
35{
36 return QStringLiteral( "rasteranalysis" );
37}
38
39void 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
66bool 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
115QVariantMap 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
230QString QgsRasterFrequencyByEqualOperatorAlgorithm::displayName() const
231{
232 return QObject::tr( "Equal to frequency" );
233}
234
235QString QgsRasterFrequencyByEqualOperatorAlgorithm::name() const
236{
237 return QStringLiteral( "equaltofrequency" );
238}
239
240QStringList QgsRasterFrequencyByEqualOperatorAlgorithm::tags() const
241{
242 return QObject::tr( "cell,equal,frequency,pixel,stack" ).split( ',' );
243}
244
245QString 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
259QgsRasterFrequencyByEqualOperatorAlgorithm *QgsRasterFrequencyByEqualOperatorAlgorithm::createInstance() const
260{
261 return new QgsRasterFrequencyByEqualOperatorAlgorithm();
262}
263
264int 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
273QString QgsRasterFrequencyByGreaterThanOperatorAlgorithm::displayName() const
274{
275 return QObject::tr( "Greater than frequency" );
276}
277
278QString QgsRasterFrequencyByGreaterThanOperatorAlgorithm::name() const
279{
280 return QStringLiteral( "greaterthanfrequency" );
281}
282
283QStringList QgsRasterFrequencyByGreaterThanOperatorAlgorithm::tags() const
284{
285 return QObject::tr( "cell,greater,frequency,pixel,stack" ).split( ',' );
286}
287
288QString 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
302QgsRasterFrequencyByGreaterThanOperatorAlgorithm *QgsRasterFrequencyByGreaterThanOperatorAlgorithm::createInstance() const
303{
304 return new QgsRasterFrequencyByGreaterThanOperatorAlgorithm();
305}
306
307int 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
316QString QgsRasterFrequencyByLessThanOperatorAlgorithm::displayName() const
317{
318 return QObject::tr( "Less than frequency" );
319}
320
321QString QgsRasterFrequencyByLessThanOperatorAlgorithm::name() const
322{
323 return QStringLiteral( "lessthanfrequency" );
324}
325
326QStringList QgsRasterFrequencyByLessThanOperatorAlgorithm::tags() const
327{
328 return QObject::tr( "cell,less,lower,frequency,pixel,stack" ).split( ',' );
329}
330
331QString 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
345QgsRasterFrequencyByLessThanOperatorAlgorithm *QgsRasterFrequencyByLessThanOperatorAlgorithm::createInstance() const
346{
347 return new QgsRasterFrequencyByLessThanOperatorAlgorithm();
348}
349
350int 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
@ 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 raster band parameter 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