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