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