QGIS API Documentation 3.43.0-Master (e01d6d7c4c0)
qgsalgorithmcellstatistics.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsalgorithmcellstatistics.cpp
3 ---------------------
4 begin : May 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
26QString QgsCellStatisticsAlgorithmBase::group() const
27{
28 return QObject::tr( "Raster analysis" );
29}
30
31QString QgsCellStatisticsAlgorithmBase::groupId() const
32{
33 return QStringLiteral( "rasteranalysis" );
34}
35
36void QgsCellStatisticsAlgorithmBase::initAlgorithm( const QVariantMap & )
37{
38 addParameter( new QgsProcessingParameterMultipleLayers( QStringLiteral( "INPUT" ), QObject::tr( "Input layers" ), Qgis::ProcessingSourceType::Raster ) );
39
40 addSpecificAlgorithmParams();
41
42 addParameter( new QgsProcessingParameterBoolean( QStringLiteral( "IGNORE_NODATA" ), QObject::tr( "Ignore NoData values" ), true ) );
43
44 addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "REFERENCE_LAYER" ), QObject::tr( "Reference layer" ) ) );
45
46 auto output_nodata_parameter = std::make_unique<QgsProcessingParameterNumber>( QStringLiteral( "OUTPUT_NODATA_VALUE" ), QObject::tr( "Output NoData value" ), Qgis::ProcessingNumberParameterType::Double, -9999, false );
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
64 addOutput( new QgsProcessingOutputString( QStringLiteral( "EXTENT" ), QObject::tr( "Extent" ) ) );
65 addOutput( new QgsProcessingOutputString( QStringLiteral( "CRS_AUTHID" ), QObject::tr( "CRS authority identifier" ) ) );
66 addOutput( new QgsProcessingOutputNumber( QStringLiteral( "WIDTH_IN_PIXELS" ), QObject::tr( "Width in pixels" ) ) );
67 addOutput( new QgsProcessingOutputNumber( QStringLiteral( "HEIGHT_IN_PIXELS" ), QObject::tr( "Height in pixels" ) ) );
68 addOutput( new QgsProcessingOutputNumber( QStringLiteral( "TOTAL_PIXEL_COUNT" ), QObject::tr( "Total pixel count" ) ) );
69}
70
71bool QgsCellStatisticsAlgorithmBase::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
72{
73 QgsRasterLayer *referenceLayer = parameterAsRasterLayer( parameters, QStringLiteral( "REFERENCE_LAYER" ), context );
74 if ( !referenceLayer )
75 throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "REFERENCE_LAYER" ) ) );
76
77 mIgnoreNoData = parameterAsBool( parameters, QStringLiteral( "IGNORE_NODATA" ), context );
78 mNoDataValue = parameterAsDouble( parameters, QStringLiteral( "OUTPUT_NODATA_VALUE" ), context );
79 mCrs = referenceLayer->crs();
80 mRasterUnitsPerPixelX = referenceLayer->rasterUnitsPerPixelX();
81 mRasterUnitsPerPixelY = referenceLayer->rasterUnitsPerPixelY();
82 mLayerWidth = referenceLayer->width();
83 mLayerHeight = referenceLayer->height();
84 mExtent = referenceLayer->extent();
85
86 const QList<QgsMapLayer *> layers = parameterAsLayerList( parameters, QStringLiteral( "INPUT" ), 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 //determine output raster data type
115 //initially raster data type to most primitive data type that is possible
116 mDataType = Qgis::DataType::Byte;
117 for ( const QgsRasterAnalysisUtils::RasterLogicInput &i : std::as_const( mInputs ) )
118 {
119 for ( int band : i.bands )
120 {
121 Qgis::DataType inputDataType = i.interface->dataType( band );
122 if ( static_cast<int>( mDataType ) < static_cast<int>( inputDataType ) )
123 mDataType = inputDataType; //if raster data type is more potent, set it as new data type
124 }
125 }
126
127 prepareSpecificAlgorithmParameters( parameters, context, feedback );
128
129 return true;
130}
131
132
133QVariantMap QgsCellStatisticsAlgorithmBase::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
134{
135 QString creationOptions = parameterAsString( parameters, QStringLiteral( "CREATION_OPTIONS" ), context ).trimmed();
136 // handle backwards compatibility parameter CREATE_OPTIONS
137 const QString optionsString = parameterAsString( parameters, QStringLiteral( "CREATE_OPTIONS" ), context );
138 if ( !optionsString.isEmpty() )
139 creationOptions = optionsString;
140
141 const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context );
142 QFileInfo fi( outputFile );
143 const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
144
145 auto writer = std::make_unique<QgsRasterFileWriter>( outputFile );
146 writer->setOutputProviderKey( QStringLiteral( "gdal" ) );
147 if ( !creationOptions.isEmpty() )
148 {
149 writer->setCreationOptions( creationOptions.split( '|' ) );
150 }
151 writer->setOutputFormat( outputFormat );
152 mOutputRasterDataProvider.reset( writer->createOneBandRaster( mDataType, mLayerWidth, mLayerHeight, mExtent, mCrs ) );
153 if ( !mOutputRasterDataProvider )
154 throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( outputFile ) );
155 if ( !mOutputRasterDataProvider->isValid() )
156 throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( outputFile, mOutputRasterDataProvider->error().message( QgsErrorMessage::Text ) ) );
157
158 mOutputRasterDataProvider->setNoDataValue( 1, mNoDataValue );
159 qgssize layerSize = static_cast<qgssize>( mLayerWidth ) * static_cast<qgssize>( mLayerHeight );
160
161 //call child statistics method
162 processRasterStack( feedback );
163
164 mOutputRasterDataProvider.reset();
165
166 QVariantMap outputs;
167 outputs.insert( QStringLiteral( "EXTENT" ), mExtent.toString() );
168 outputs.insert( QStringLiteral( "CRS_AUTHID" ), mCrs.authid() );
169 outputs.insert( QStringLiteral( "WIDTH_IN_PIXELS" ), mLayerWidth );
170 outputs.insert( QStringLiteral( "HEIGHT_IN_PIXELS" ), mLayerHeight );
171 outputs.insert( QStringLiteral( "TOTAL_PIXEL_COUNT" ), layerSize );
172 outputs.insert( QStringLiteral( "OUTPUT" ), outputFile );
173
174 return outputs;
175}
176
177
178//
179//QgsCellStatisticsAlgorithm
180//
181QString QgsCellStatisticsAlgorithm::displayName() const
182{
183 return QObject::tr( "Cell statistics" );
184}
185
186QString QgsCellStatisticsAlgorithm::name() const
187{
188 return QStringLiteral( "cellstatistics" );
189}
190
191QStringList QgsCellStatisticsAlgorithm::tags() const
192{
193 return QObject::tr( "cell,pixel,statistic,count,mean,sum,majority,minority,variance,variety,range,median,minimum,maximum" ).split( ',' );
194}
195
196QString QgsCellStatisticsAlgorithm::shortHelpString() const
197{
198 return QObject::tr( "The Cell statistics algorithm computes a value for each cell of the "
199 "output raster. At each cell location, "
200 "the output value is defined as a function of all overlaid cell values of the "
201 "input rasters.\n\n"
202 "The output raster's extent and resolution is defined by a reference "
203 "raster. The following functions can be applied on the input "
204 "raster cells per output raster cell location:\n"
205 "<ul> "
206 " <li>Sum</li>"
207 " <li>Count</li>"
208 " <li>Mean</li>"
209 " <li>Median</li>"
210 " <li>Standard deviation</li>"
211 " <li>Variance</li>"
212 " <li>Minimum</li>"
213 " <li>Maximum</li>"
214 " <li>Minority (least frequent value)</li>"
215 " <li>Majority (most frequent value)</li>"
216 " <li>Range (max-min)</li>"
217 " <li>Variety (count of unique values)</li>"
218 "</ul> "
219 "Input raster layers that do not match the cell size of the reference raster layer will be "
220 "resampled using nearest neighbor resampling. The output raster data type will be set to "
221 "the most complex data type present in the input datasets except when using the functions "
222 "Mean, Standard deviation and Variance (data type is always Float32/Float64 depending on input float type) or Count and Variety (data type is always Int32).\n"
223 "<i>Calculation details - general:</i> NoData values in any of the input layers will result in a NoData cell output if the Ignore NoData parameter is not set.\n"
224 "<i>Calculation details - Count:</i> Count will always result in the number of cells without NoData values at the current cell location.\n"
225 "<i>Calculation details - Median:</i> If the number of input layers is even, the median will be calculated as the "
226 "arithmetic mean of the two middle values of the ordered cell input values. In this case the output data type is Float32.\n"
227 "<i>Calculation details - Minority/Majority:</i> If no unique minority or majority could be found, the result is NoData, except all "
228 "input cell values are equal." );
229}
230
231QString QgsCellStatisticsAlgorithm::shortDescription() const
232{
233 return QObject::tr( "Generates a raster whose cell values are computed from overlaid cell values of the input rasters." );
234}
235
236QgsCellStatisticsAlgorithm *QgsCellStatisticsAlgorithm::createInstance() const
237{
238 return new QgsCellStatisticsAlgorithm();
239}
240
241void QgsCellStatisticsAlgorithm::addSpecificAlgorithmParams()
242{
243 QStringList statistics = QStringList();
244 statistics << QObject::tr( "Sum" )
245 << QObject::tr( "Count" )
246 << QObject::tr( "Mean" )
247 << QObject::tr( "Median" )
248 << QObject::tr( "Standard deviation" )
249 << QObject::tr( "Variance" )
250 << QObject::tr( "Minimum" )
251 << QObject::tr( "Maximum" )
252 << QObject::tr( "Minority" )
253 << QObject::tr( "Majority" )
254 << QObject::tr( "Range" )
255 << QObject::tr( "Variety" );
256
257 addParameter( new QgsProcessingParameterEnum( QStringLiteral( "STATISTIC" ), QObject::tr( "Statistic" ), statistics, false, 0, false ) );
258}
259
260bool QgsCellStatisticsAlgorithm::prepareSpecificAlgorithmParameters( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
261{
262 Q_UNUSED( feedback )
263 //obtain statistic method
264 mMethod = static_cast<QgsRasterAnalysisUtils::CellValueStatisticMethods>( parameterAsEnum( parameters, QStringLiteral( "STATISTIC" ), context ) );
265
266 //force data types on specific functions in the cellstatistics alg if input data types don't match
267 if (
268 mMethod == QgsRasterAnalysisUtils::Mean || mMethod == QgsRasterAnalysisUtils::StandardDeviation || mMethod == QgsRasterAnalysisUtils::Variance || ( mMethod == QgsRasterAnalysisUtils::Median && ( mInputs.size() % 2 == 0 ) )
269 )
270 {
271 if ( static_cast<int>( mDataType ) < 6 )
272 mDataType = Qgis::DataType::Float32; //force float on mean, stddev and median with equal number of input layers if all inputs are integer
273 }
274 else if ( mMethod == QgsRasterAnalysisUtils::Count || mMethod == QgsRasterAnalysisUtils::Variety ) //count, variety
275 {
276 if ( static_cast<int>( mDataType ) > 5 ) //if is floating point type
277 mDataType = Qgis::DataType::Int32; //force integer on variety if all inputs are float or complex
278 }
279 return true;
280}
281
282void QgsCellStatisticsAlgorithm::processRasterStack( QgsProcessingFeedback *feedback )
283{
286 int nbBlocksWidth = static_cast<int>( std::ceil( 1.0 * mLayerWidth / maxWidth ) );
287 int nbBlocksHeight = static_cast<int>( std::ceil( 1.0 * mLayerHeight / maxHeight ) );
288 int nbBlocks = nbBlocksWidth * nbBlocksHeight;
289 mOutputRasterDataProvider->setEditable( true );
290 QgsRasterIterator outputIter( mOutputRasterDataProvider.get() );
291 outputIter.startRasterRead( 1, mLayerWidth, mLayerHeight, mExtent );
292
293 int iterLeft = 0;
294 int iterTop = 0;
295 int iterCols = 0;
296 int iterRows = 0;
297 QgsRectangle blockExtent;
298 std::unique_ptr<QgsRasterBlock> outputBlock;
299 while ( outputIter.readNextRasterPart( 1, iterCols, iterRows, outputBlock, iterLeft, iterTop, &blockExtent ) )
300 {
301 std::vector<std::unique_ptr<QgsRasterBlock>> inputBlocks;
302 for ( const QgsRasterAnalysisUtils::RasterLogicInput &i : std::as_const( mInputs ) )
303 {
304 if ( feedback->isCanceled() )
305 break; //in case some slow data sources are loaded
306 for ( int band : i.bands )
307 {
308 if ( feedback->isCanceled() )
309 break; //in case some slow data sources are loaded
310 std::unique_ptr<QgsRasterBlock> b( i.interface->block( band, blockExtent, iterCols, iterRows ) );
311 inputBlocks.emplace_back( std::move( b ) );
312 }
313 }
314
315 feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
316 for ( int row = 0; row < iterRows; row++ )
317 {
318 if ( feedback->isCanceled() )
319 break;
320
321 for ( int col = 0; col < iterCols; col++ )
322 {
323 double result = 0;
324 bool noDataInStack = false;
325 std::vector<double> cellValues = QgsRasterAnalysisUtils::getCellValuesFromBlockStack( inputBlocks, row, col, noDataInStack );
326 int cellValueStackSize = cellValues.size();
327
328 if ( noDataInStack && !mIgnoreNoData )
329 {
330 //output cell will always be NoData if NoData occurs in cellValueStack and NoData is not ignored
331 //this saves unnecessary iterations on the cellValueStack
332 if ( mMethod == QgsRasterAnalysisUtils::Count )
333 outputBlock->setValue( row, col, cellValueStackSize );
334 else
335 {
336 outputBlock->setValue( row, col, mNoDataValue );
337 }
338 }
339 else if ( !noDataInStack || ( mIgnoreNoData && cellValueStackSize > 0 ) )
340 {
341 switch ( mMethod )
342 {
343 case QgsRasterAnalysisUtils::Sum:
344 result = std::accumulate( cellValues.begin(), cellValues.end(), 0.0 );
345 break;
346 case QgsRasterAnalysisUtils::Count:
347 result = cellValueStackSize;
348 break;
349 case QgsRasterAnalysisUtils::Mean:
350 result = QgsRasterAnalysisUtils::meanFromCellValues( cellValues, cellValueStackSize );
351 break;
352 case QgsRasterAnalysisUtils::Median:
353 result = QgsRasterAnalysisUtils::medianFromCellValues( cellValues, cellValueStackSize );
354 break;
355 case QgsRasterAnalysisUtils::StandardDeviation:
356 result = QgsRasterAnalysisUtils::stddevFromCellValues( cellValues, cellValueStackSize );
357 break;
358 case QgsRasterAnalysisUtils::Variance:
359 result = QgsRasterAnalysisUtils::varianceFromCellValues( cellValues, cellValueStackSize );
360 break;
361 case QgsRasterAnalysisUtils::Minimum:
362 result = QgsRasterAnalysisUtils::minimumFromCellValues( cellValues );
363 break;
364 case QgsRasterAnalysisUtils::Maximum:
365 result = QgsRasterAnalysisUtils::maximumFromCellValues( cellValues );
366 break;
367 case QgsRasterAnalysisUtils::Minority:
368 result = QgsRasterAnalysisUtils::minorityFromCellValues( cellValues, mNoDataValue, cellValueStackSize );
369 break;
370 case QgsRasterAnalysisUtils::Majority:
371 result = QgsRasterAnalysisUtils::majorityFromCellValues( cellValues, mNoDataValue, cellValueStackSize );
372 break;
373 case QgsRasterAnalysisUtils::Range:
374 result = QgsRasterAnalysisUtils::rangeFromCellValues( cellValues );
375 break;
376 case QgsRasterAnalysisUtils::Variety:
377 result = QgsRasterAnalysisUtils::varietyFromCellValues( cellValues );
378 break;
379 }
380 outputBlock->setValue( row, col, result );
381 }
382 else
383 {
384 //result is NoData if cellValueStack contains no valid values, eg. all cellValues are NoData
385 outputBlock->setValue( row, col, mNoDataValue );
386 }
387 }
388 }
389 if ( !mOutputRasterDataProvider->writeBlock( outputBlock.get(), 1, iterLeft, iterTop ) )
390 {
391 throw QgsProcessingException( QObject::tr( "Could not write raster block: %1" ).arg( mOutputRasterDataProvider->error().summary() ) );
392 }
393 }
394 mOutputRasterDataProvider->setEditable( false );
395}
396
397//
398//QgsCellStatisticsPercentileAlgorithm
399//
400QString QgsCellStatisticsPercentileAlgorithm::displayName() const
401{
402 return QObject::tr( "Cell stack percentile" );
403}
404
405QString QgsCellStatisticsPercentileAlgorithm::name() const
406{
407 return QStringLiteral( "cellstackpercentile" );
408}
409
410QStringList QgsCellStatisticsPercentileAlgorithm::tags() const
411{
412 return QObject::tr( "cell,pixel,statistic,percentile,quantile,quartile" ).split( ',' );
413}
414
415QString QgsCellStatisticsPercentileAlgorithm::shortHelpString() const
416{
417 return QObject::tr( "This algorithm generates a raster containing the cell-wise percentile value of a stack of input rasters. "
418 "The percentile to return is determined by the percentile input value (ranges between 0 and 1). "
419 "At each cell location, the specified percentile is obtained using the respective value from "
420 "the stack of all overlaid and sorted cell values of the input rasters.\n\n"
421 "There are three methods for percentile calculation:"
422 "<ul> "
423 " <li>Nearest rank</li>"
424 " <li>Inclusive linear interpolation (PERCENTILE.INC)</li>"
425 " <li>Exclusive linear interpolation (PERCENTILE.EXC)</li>"
426 "</ul> "
427 "While the output value can stay the same for the nearest rank method (obtains the value that is nearest to the "
428 "specified percentile), the linear interpolation method return unique values for different percentiles. Both interpolation "
429 "methods follow their counterpart methods implemented by LibreOffice or Microsoft Excel. \n\n"
430 "The output raster's extent and resolution is defined by a reference "
431 "raster. If the input raster layers that do not match the cell size of the reference raster layer will be "
432 "resampled using nearest neighbor resampling. NoData values in any of the input layers will result in a NoData cell output if the Ignore NoData parameter is not set. "
433 "The output raster data type will be set to the most complex data type present in the input datasets. " );
434}
435
436QString QgsCellStatisticsPercentileAlgorithm::shortDescription() const
437{
438 return QObject::tr( "Generates a raster containing the cell-wise percentile value of a stack of input rasters." );
439}
440
441QgsCellStatisticsPercentileAlgorithm *QgsCellStatisticsPercentileAlgorithm::createInstance() const
442{
443 return new QgsCellStatisticsPercentileAlgorithm();
444}
445
446void QgsCellStatisticsPercentileAlgorithm::addSpecificAlgorithmParams()
447{
448 addParameter( new QgsProcessingParameterEnum( QStringLiteral( "METHOD" ), QObject::tr( "Method" ), QStringList() << QObject::tr( "Nearest rank" ) << QObject::tr( "Inclusive linear interpolation (PERCENTILE.INC)" ) << QObject::tr( "Exclusive linear interpolation (PERCENTILE.EXC)" ), false, 0, false ) );
449 addParameter( new QgsProcessingParameterNumber( QStringLiteral( "PERCENTILE" ), QObject::tr( "Percentile" ), Qgis::ProcessingNumberParameterType::Double, 0.25, false, 0.0, 1.0 ) );
450}
451
452bool QgsCellStatisticsPercentileAlgorithm::prepareSpecificAlgorithmParameters( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
453{
454 Q_UNUSED( feedback )
455 mMethod = static_cast<QgsRasterAnalysisUtils::CellValuePercentileMethods>( parameterAsEnum( parameters, QStringLiteral( "METHOD" ), context ) );
456 mPercentile = parameterAsDouble( parameters, QStringLiteral( "PERCENTILE" ), context );
457
458 //default percentile output data type to float32 raster if interpolation method is chosen
459 //otherwise use the most potent data type in the input raster stack (see prepareAlgorithm() in base class)
460 if ( mMethod != QgsRasterAnalysisUtils::CellValuePercentileMethods::NearestRankPercentile && static_cast<int>( mDataType ) < 6 )
461 mDataType = Qgis::DataType::Float32;
462
463 return true;
464}
465
466void QgsCellStatisticsPercentileAlgorithm::processRasterStack( QgsProcessingFeedback *feedback )
467{
470 int nbBlocksWidth = static_cast<int>( std::ceil( 1.0 * mLayerWidth / maxWidth ) );
471 int nbBlocksHeight = static_cast<int>( std::ceil( 1.0 * mLayerHeight / maxHeight ) );
472 int nbBlocks = nbBlocksWidth * nbBlocksHeight;
473 mOutputRasterDataProvider->setEditable( true );
474 QgsRasterIterator outputIter( mOutputRasterDataProvider.get() );
475 outputIter.startRasterRead( 1, mLayerWidth, mLayerHeight, mExtent );
476
477 int iterLeft = 0;
478 int iterTop = 0;
479 int iterCols = 0;
480 int iterRows = 0;
481 QgsRectangle blockExtent;
482 std::unique_ptr<QgsRasterBlock> outputBlock;
483 while ( outputIter.readNextRasterPart( 1, iterCols, iterRows, outputBlock, iterLeft, iterTop, &blockExtent ) )
484 {
485 std::vector<std::unique_ptr<QgsRasterBlock>> inputBlocks;
486 for ( const QgsRasterAnalysisUtils::RasterLogicInput &i : std::as_const( mInputs ) )
487 {
488 if ( feedback->isCanceled() )
489 break; //in case some slow data sources are loaded
490 for ( int band : i.bands )
491 {
492 if ( feedback->isCanceled() )
493 break; //in case some slow data sources are loaded
494 std::unique_ptr<QgsRasterBlock> b( i.interface->block( band, blockExtent, iterCols, iterRows ) );
495 inputBlocks.emplace_back( std::move( b ) );
496 }
497 }
498
499 feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
500 for ( int row = 0; row < iterRows; row++ )
501 {
502 if ( feedback->isCanceled() )
503 break;
504
505 for ( int col = 0; col < iterCols; col++ )
506 {
507 double result = 0;
508 bool noDataInStack = false;
509 std::vector<double> cellValues = QgsRasterAnalysisUtils::getCellValuesFromBlockStack( inputBlocks, row, col, noDataInStack );
510 int cellValueStackSize = cellValues.size();
511
512 if ( noDataInStack && !mIgnoreNoData )
513 {
514 outputBlock->setValue( row, col, mNoDataValue );
515 }
516 else if ( !noDataInStack || ( mIgnoreNoData && cellValueStackSize > 0 ) )
517 {
518 switch ( mMethod )
519 {
520 case QgsRasterAnalysisUtils::NearestRankPercentile:
521 result = QgsRasterAnalysisUtils::nearestRankPercentile( cellValues, cellValueStackSize, mPercentile );
522 break;
523 case QgsRasterAnalysisUtils::InterpolatedPercentileInc:
524 result = QgsRasterAnalysisUtils::interpolatedPercentileInc( cellValues, cellValueStackSize, mPercentile );
525 break;
526 case QgsRasterAnalysisUtils::InterpolatedPercentileExc:
527 result = QgsRasterAnalysisUtils::interpolatedPercentileExc( cellValues, cellValueStackSize, mPercentile, mNoDataValue );
528 break;
529 }
530 outputBlock->setValue( row, col, result );
531 }
532 else
533 {
534 //result is NoData if cellValueStack contains no valid values, eg. all cellValues are NoData
535 outputBlock->setValue( row, col, mNoDataValue );
536 }
537 }
538 }
539 if ( !mOutputRasterDataProvider->writeBlock( outputBlock.get(), 1, iterLeft, iterTop ) )
540 {
541 throw QgsProcessingException( QObject::tr( "Could not write raster block: %1" ).arg( mOutputRasterDataProvider->error().summary() ) );
542 }
543 }
544 mOutputRasterDataProvider->setEditable( false );
545}
546
547//
548//QgsCellStatisticsPercentRankFromValueAlgorithm
549//
550QString QgsCellStatisticsPercentRankFromValueAlgorithm::displayName() const
551{
552 return QObject::tr( "Cell stack percent rank from value" );
553}
554
555QString QgsCellStatisticsPercentRankFromValueAlgorithm::name() const
556{
557 return QStringLiteral( "cellstackpercentrankfromvalue" );
558}
559
560QStringList QgsCellStatisticsPercentRankFromValueAlgorithm::tags() const
561{
562 return QObject::tr( "cell,pixel,statistic,percentrank,rank,percent,value" ).split( ',' );
563}
564
565QString QgsCellStatisticsPercentRankFromValueAlgorithm::shortHelpString() const
566{
567 return QObject::tr( "This algorithm generates a raster containing the cell-wise percent rank value of a stack of input rasters based on a single input value.\n\n"
568 "At each cell location, the specified value is ranked among the respective values in the stack of all overlaid and sorted cell values from the input rasters. "
569 "For values outside of the stack value distribution, the algorithm returns NoData because the value cannot be ranked among the cell values.\n\n"
570 "There are two methods for percentile calculation:"
571 "<ul> "
572 " <li>Inclusive linearly interpolated percent rank (PERCENTRANK.INC)</li>"
573 " <li>Exclusive linearly interpolated percent rank (PERCENTRANK.EXC)</li>"
574 "</ul> "
575 "The linear interpolation method return the unique percent rank for different values. Both interpolation "
576 "methods follow their counterpart methods implemented by LibreOffice or Microsoft Excel. \n\n"
577 "The output raster's extent and resolution is defined by a reference "
578 "raster. If the input raster layers that do not match the cell size of the reference raster layer will be "
579 "resampled using nearest neighbor resampling. NoData values in any of the input layers will result in a NoData cell output if the Ignore NoData parameter is not set. "
580 "The output raster data type will always be Float32." );
581}
582
583QString QgsCellStatisticsPercentRankFromValueAlgorithm::shortDescription() const
584{
585 return QObject::tr( "Generates a raster containing the cell-wise percent rank value of a stack of input rasters based on a single input value." );
586}
587
588QgsCellStatisticsPercentRankFromValueAlgorithm *QgsCellStatisticsPercentRankFromValueAlgorithm::createInstance() const
589{
590 return new QgsCellStatisticsPercentRankFromValueAlgorithm();
591}
592
593void QgsCellStatisticsPercentRankFromValueAlgorithm::addSpecificAlgorithmParams()
594{
595 addParameter( new QgsProcessingParameterEnum( QStringLiteral( "METHOD" ), QObject::tr( "Method" ), QStringList() << QObject::tr( "Inclusive linear interpolation (PERCENTRANK.INC)" ) << QObject::tr( "Exclusive linear interpolation (PERCENTRANK.EXC)" ), false, 0, false ) );
596 addParameter( new QgsProcessingParameterNumber( QStringLiteral( "VALUE" ), QObject::tr( "Value" ), Qgis::ProcessingNumberParameterType::Double, 10, false ) );
597}
598
599bool QgsCellStatisticsPercentRankFromValueAlgorithm::prepareSpecificAlgorithmParameters( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
600{
601 Q_UNUSED( feedback )
602 mMethod = static_cast<QgsRasterAnalysisUtils::CellValuePercentRankMethods>( parameterAsEnum( parameters, QStringLiteral( "METHOD" ), context ) );
603 mValue = parameterAsDouble( parameters, QStringLiteral( "VALUE" ), context );
604
605 //output data type always defaults to Float32 because result only ranges between 0 and 1
606 mDataType = Qgis::DataType::Float32;
607 return true;
608}
609
610void QgsCellStatisticsPercentRankFromValueAlgorithm::processRasterStack( QgsProcessingFeedback *feedback )
611{
614 int nbBlocksWidth = static_cast<int>( std::ceil( 1.0 * mLayerWidth / maxWidth ) );
615 int nbBlocksHeight = static_cast<int>( std::ceil( 1.0 * mLayerHeight / maxHeight ) );
616 int nbBlocks = nbBlocksWidth * nbBlocksHeight;
617 mOutputRasterDataProvider->setEditable( true );
618 QgsRasterIterator outputIter( mOutputRasterDataProvider.get() );
619 outputIter.startRasterRead( 1, mLayerWidth, mLayerHeight, mExtent );
620
621 int iterLeft = 0;
622 int iterTop = 0;
623 int iterCols = 0;
624 int iterRows = 0;
625 QgsRectangle blockExtent;
626 std::unique_ptr<QgsRasterBlock> outputBlock;
627 while ( outputIter.readNextRasterPart( 1, iterCols, iterRows, outputBlock, iterLeft, iterTop, &blockExtent ) )
628 {
629 std::vector<std::unique_ptr<QgsRasterBlock>> inputBlocks;
630 for ( const QgsRasterAnalysisUtils::RasterLogicInput &i : std::as_const( mInputs ) )
631 {
632 if ( feedback->isCanceled() )
633 break; //in case some slow data sources are loaded
634 for ( int band : i.bands )
635 {
636 if ( feedback->isCanceled() )
637 break; //in case some slow data sources are loaded
638 std::unique_ptr<QgsRasterBlock> b( i.interface->block( band, blockExtent, iterCols, iterRows ) );
639 inputBlocks.emplace_back( std::move( b ) );
640 }
641 }
642
643 feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
644 for ( int row = 0; row < iterRows; row++ )
645 {
646 if ( feedback->isCanceled() )
647 break;
648
649 for ( int col = 0; col < iterCols; col++ )
650 {
651 double result = 0;
652 bool noDataInStack = false;
653 std::vector<double> cellValues = QgsRasterAnalysisUtils::getCellValuesFromBlockStack( inputBlocks, row, col, noDataInStack );
654 int cellValueStackSize = cellValues.size();
655
656 if ( noDataInStack && !mIgnoreNoData )
657 {
658 outputBlock->setValue( row, col, mNoDataValue );
659 }
660 else if ( !noDataInStack || ( mIgnoreNoData && cellValueStackSize > 0 ) )
661 {
662 switch ( mMethod )
663 {
664 case QgsRasterAnalysisUtils::InterpolatedPercentRankInc:
665 result = QgsRasterAnalysisUtils::interpolatedPercentRankInc( cellValues, cellValueStackSize, mValue, mNoDataValue );
666 break;
667 case QgsRasterAnalysisUtils::InterpolatedPercentRankExc:
668 result = QgsRasterAnalysisUtils::interpolatedPercentRankExc( cellValues, cellValueStackSize, mValue, mNoDataValue );
669 break;
670 }
671 outputBlock->setValue( row, col, result );
672 }
673 else
674 {
675 //result is NoData if cellValueStack contains no valid values, eg. all cellValues are NoData
676 outputBlock->setValue( row, col, mNoDataValue );
677 }
678 }
679 }
680 if ( !mOutputRasterDataProvider->writeBlock( outputBlock.get(), 1, iterLeft, iterTop ) )
681 {
682 throw QgsProcessingException( QObject::tr( "Could not write raster block: %1" ).arg( mOutputRasterDataProvider->error().summary() ) );
683 }
684 }
685 mOutputRasterDataProvider->setEditable( false );
686}
687
688
689//
690//QgsCellStatisticsPercentRankFromRasterAlgorithm
691//
692QString QgsCellStatisticsPercentRankFromRasterAlgorithm::displayName() const
693{
694 return QObject::tr( "Cell stack percentrank from raster layer" );
695}
696
697QString QgsCellStatisticsPercentRankFromRasterAlgorithm::name() const
698{
699 return QStringLiteral( "cellstackpercentrankfromrasterlayer" );
700}
701
702QStringList QgsCellStatisticsPercentRankFromRasterAlgorithm::tags() const
703{
704 return QObject::tr( "cell,pixel,statistic,percentrank,rank,percent,value,raster" ).split( ',' );
705}
706
707QString QgsCellStatisticsPercentRankFromRasterAlgorithm::shortHelpString() const
708{
709 return QObject::tr( "This algorithm generates a raster containing the cell-wise percent rank value of a stack of input rasters "
710 "based on an input value raster.\n\n"
711 "At each cell location, the current value of the value raster is used ranked among the respective values in the stack of all overlaid and sorted cell values of the input rasters. "
712 "For values outside of the the stack value distribution, the algorithm returns NoData because the value cannot be ranked among the cell values.\n\n"
713 "There are two methods for percentile calculation:"
714 "<ul> "
715 " <li>Inclusive linearly interpolated percent rank (PERCENTRANK.INC)</li>"
716 " <li>Exclusive linearly interpolated percent rank (PERCENTRANK.EXC)</li>"
717 "</ul> "
718 "The linear interpolation method return the unique percent rank for different values. Both interpolation "
719 "methods follow their counterpart methods implemented by LibreOffice or Microsoft Excel. \n\n"
720 "The output raster's extent and resolution is defined by a reference "
721 "raster. If the input raster layers that do not match the cell size of the reference raster layer will be "
722 "resampled using nearest neighbor resampling. NoData values in any of the input layers will result in a NoData cell output if the Ignore NoData parameter is not set. "
723 "The output raster data type will always be Float32." );
724}
725
726QString QgsCellStatisticsPercentRankFromRasterAlgorithm::shortDescription() const
727{
728 return QObject::tr( "Generates a raster containing the cell-wise percent rank value of a stack of input rasters based on an input value raster." );
729}
730
731QgsCellStatisticsPercentRankFromRasterAlgorithm *QgsCellStatisticsPercentRankFromRasterAlgorithm::createInstance() const
732{
733 return new QgsCellStatisticsPercentRankFromRasterAlgorithm();
734}
735
736void QgsCellStatisticsPercentRankFromRasterAlgorithm::addSpecificAlgorithmParams()
737{
738 addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT_VALUE_RASTER" ), QObject::tr( "Value raster layer" ) ) );
739 addParameter( new QgsProcessingParameterBand( QStringLiteral( "VALUE_RASTER_BAND" ), QObject::tr( "Value raster band" ), 1, QStringLiteral( "VALUE_LAYER" ) ) );
740 addParameter( new QgsProcessingParameterEnum( QStringLiteral( "METHOD" ), QObject::tr( "Method" ), QStringList() << QObject::tr( "Inclusive linear interpolation (PERCENTRANK.INC)" ) << QObject::tr( "Exclusive linear interpolation (PERCENTRANK.EXC)" ), false, 0, false ) );
741}
742
743bool QgsCellStatisticsPercentRankFromRasterAlgorithm::prepareSpecificAlgorithmParameters( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
744{
745 Q_UNUSED( feedback )
746 mMethod = static_cast<QgsRasterAnalysisUtils::CellValuePercentRankMethods>( parameterAsEnum( parameters, QStringLiteral( "METHOD" ), context ) );
747
748 QgsRasterLayer *inputValueRaster = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT_VALUE_RASTER" ), context );
749 if ( !inputValueRaster )
750 throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT_VALUE_RASTER" ) ) );
751
752 mValueRasterInterface.reset( inputValueRaster->dataProvider()->clone() );
753
754 mValueRasterBand = parameterAsInt( parameters, QStringLiteral( "VALUE_RASTER_BAND" ), context );
755
756 //output data type always defaults to Float32 because result only ranges between 0 and 1
757 mDataType = Qgis::DataType::Float32;
758 return true;
759}
760
761void QgsCellStatisticsPercentRankFromRasterAlgorithm::processRasterStack( QgsProcessingFeedback *feedback )
762{
765 int nbBlocksWidth = static_cast<int>( std::ceil( 1.0 * mLayerWidth / maxWidth ) );
766 int nbBlocksHeight = static_cast<int>( std::ceil( 1.0 * mLayerHeight / maxHeight ) );
767 int nbBlocks = nbBlocksWidth * nbBlocksHeight;
768 mOutputRasterDataProvider->setEditable( true );
769 QgsRasterIterator outputIter( mOutputRasterDataProvider.get() );
770 outputIter.startRasterRead( 1, mLayerWidth, mLayerHeight, mExtent );
771
772 int iterLeft = 0;
773 int iterTop = 0;
774 int iterCols = 0;
775 int iterRows = 0;
776 QgsRectangle blockExtent;
777 std::unique_ptr<QgsRasterBlock> outputBlock;
778 while ( outputIter.readNextRasterPart( 1, iterCols, iterRows, outputBlock, iterLeft, iterTop, &blockExtent ) )
779 {
780 std::unique_ptr<QgsRasterBlock> valueBlock( mValueRasterInterface->block( mValueRasterBand, blockExtent, iterCols, iterRows ) );
781
782 std::vector<std::unique_ptr<QgsRasterBlock>> inputBlocks;
783 for ( const QgsRasterAnalysisUtils::RasterLogicInput &i : std::as_const( mInputs ) )
784 {
785 if ( feedback->isCanceled() )
786 break; //in case some slow data sources are loaded
787 for ( int band : i.bands )
788 {
789 if ( feedback->isCanceled() )
790 break; //in case some slow data sources are loaded
791 std::unique_ptr<QgsRasterBlock> b( i.interface->block( band, blockExtent, iterCols, iterRows ) );
792 inputBlocks.emplace_back( std::move( b ) );
793 }
794 }
795
796 feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
797 for ( int row = 0; row < iterRows; row++ )
798 {
799 if ( feedback->isCanceled() )
800 break;
801
802 for ( int col = 0; col < iterCols; col++ )
803 {
804 bool percentRankValueIsNoData = false;
805 double percentRankValue = valueBlock->valueAndNoData( row, col, percentRankValueIsNoData );
806
807 double result = 0;
808 bool noDataInStack = false;
809 std::vector<double> cellValues = QgsRasterAnalysisUtils::getCellValuesFromBlockStack( inputBlocks, row, col, noDataInStack );
810 int cellValueStackSize = cellValues.size();
811
812 if ( noDataInStack && !mIgnoreNoData && !percentRankValueIsNoData )
813 {
814 outputBlock->setValue( row, col, mNoDataValue );
815 }
816 else if ( !noDataInStack || ( !percentRankValueIsNoData && mIgnoreNoData && cellValueStackSize > 0 ) )
817 {
818 switch ( mMethod )
819 {
820 case QgsRasterAnalysisUtils::InterpolatedPercentRankInc:
821 result = QgsRasterAnalysisUtils::interpolatedPercentRankInc( cellValues, cellValueStackSize, percentRankValue, mNoDataValue );
822 break;
823 case QgsRasterAnalysisUtils::InterpolatedPercentRankExc:
824 result = QgsRasterAnalysisUtils::interpolatedPercentRankExc( cellValues, cellValueStackSize, percentRankValue, mNoDataValue );
825 break;
826 }
827 outputBlock->setValue( row, col, result );
828 }
829 else
830 {
831 //result is NoData if cellValueStack contains no valid values, eg. all cellValues are NoData or percentRankValue is NoData
832 outputBlock->setValue( row, col, mNoDataValue );
833 }
834 }
835 }
836 if ( !mOutputRasterDataProvider->writeBlock( outputBlock.get(), 1, iterLeft, iterTop ) )
837 {
838 throw QgsProcessingException( QObject::tr( "Could not write raster block: %1" ).arg( mOutputRasterDataProvider->error().summary() ) );
839 }
840 }
841 mOutputRasterDataProvider->setEditable( false );
842}
843
DataType
Raster data types.
Definition qgis.h:351
@ Float32
Thirty two bit floating point (float)
@ Byte
Eight bit unsigned integer (quint8)
@ 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.
An enum based parameter for processing algorithms, allowing for selection from predefined values.
A parameter for processing algorithms which accepts multiple map layers.
A numeric parameter for processing algorithms.
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:6826