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