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