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