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