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