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