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