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