QGIS API Documentation  3.23.0-Master (7c4a6de034)
qgsalgorithmroundrastervalues.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsalgorithmroundrastervalues.cpp
3  ---------------------
4  begin : April 2020
5  copyright : (C) 2020 by Clemens Raffler
6  email : clemens dot raffler at gmail dot com
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
19 #include "qgsrasterfilewriter.h"
20 
22 
23 QString QgsRoundRasterValuesAlgorithm::name() const
24 {
25  return QStringLiteral( "roundrastervalues" );
26 }
27 
28 QString QgsRoundRasterValuesAlgorithm::displayName() const
29 {
30  return QObject::tr( "Round raster" );
31 }
32 
33 QStringList QgsRoundRasterValuesAlgorithm::tags() const
34 {
35  return QObject::tr( "data,cells,round,truncate" ).split( ',' );
36 }
37 
38 QString QgsRoundRasterValuesAlgorithm::group() const
39 {
40  return QObject::tr( "Raster analysis" );
41 }
42 
43 QString QgsRoundRasterValuesAlgorithm::groupId() const
44 {
45  return QStringLiteral( "rasteranalysis" );
46 }
47 
48 void QgsRoundRasterValuesAlgorithm::initAlgorithm( const QVariantMap & )
49 {
50  addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT" ), QStringLiteral( "Input raster" ) ) );
51  addParameter( new QgsProcessingParameterBand( QStringLiteral( "BAND" ), QObject::tr( "Band number" ), 1, QStringLiteral( "INPUT" ) ) );
52  addParameter( new QgsProcessingParameterEnum( QStringLiteral( "ROUNDING_DIRECTION" ), QObject::tr( "Rounding direction" ), QStringList() << QObject::tr( "Round up" ) << QObject::tr( "Round to nearest" ) << QObject::tr( "Round down" ), false, 1 ) );
53  addParameter( new QgsProcessingParameterNumber( QStringLiteral( "DECIMAL_PLACES" ), QObject::tr( "Number of decimals places" ), QgsProcessingParameterNumber::Integer, 2 ) );
54  addParameter( new QgsProcessingParameterRasterDestination( QStringLiteral( "OUTPUT" ), QObject::tr( "Output raster" ) ) );
55  std::unique_ptr< QgsProcessingParameterDefinition > baseParameter = std::make_unique< QgsProcessingParameterNumber >( QStringLiteral( "BASE_N" ), QObject::tr( "Base n for rounding to multiples of n" ), QgsProcessingParameterNumber::Integer, 10, true, 1 );
56  baseParameter->setFlags( QgsProcessingParameterDefinition::FlagAdvanced );
57  addParameter( baseParameter.release() );
58 }
59 
60 QString QgsRoundRasterValuesAlgorithm::shortHelpString() const
61 {
62  return QObject::tr( "This algorithm rounds the cell values of a raster dataset according to the specified number of decimals.\n "
63  "Alternatively, a negative number of decimal places may be used to round values to powers of a base n "
64  "(specified in the advanced parameter Base n). For example, with a Base value n of 10 and Decimal places of -1 "
65  "the algorithm rounds cell values to multiples of 10, -2 rounds to multiples of 100, and so on. Arbitrary base values "
66  "may be chosen, the algorithm applies the same multiplicative principle. Rounding cell values to multiples of "
67  "a base n may be used to generalize raster layers.\n"
68  "The algorithm preserves the data type of the input raster. Therefore byte/integer rasters can only be rounded "
69  "to multiples of a base n, otherwise a warning is raised and the raster gets copied as byte/integer raster" );
70 }
71 
72 QgsRoundRasterValuesAlgorithm *QgsRoundRasterValuesAlgorithm::createInstance() const
73 {
74  return new QgsRoundRasterValuesAlgorithm();
75 }
76 
77 bool QgsRoundRasterValuesAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
78 {
79  Q_UNUSED( feedback );
80  QgsRasterLayer *inputRaster = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT" ), context );
81  mDecimalPrecision = parameterAsInt( parameters, QStringLiteral( "DECIMAL_PLACES" ), context );
82  mBaseN = parameterAsInt( parameters, QStringLiteral( "BASE_N" ), context );
83  mMultipleOfBaseN = pow( mBaseN, abs( mDecimalPrecision ) );
84  mScaleFactor = std::pow( 10.0, mDecimalPrecision );
85 
86  if ( !inputRaster )
87  throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT" ) ) );
88 
89  mBand = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
90  if ( mBand < 1 || mBand > inputRaster->bandCount() )
91  throw QgsProcessingException( QObject::tr( "Invalid band number for BAND (%1): Valid values for input raster are 1 to %2" ).arg( mBand ).arg( inputRaster->bandCount() ) );
92 
93  mRoundingDirection = parameterAsEnum( parameters, QStringLiteral( "ROUNDING_DIRECTION" ), context );
94 
95  mInterface.reset( inputRaster->dataProvider()->clone() );
96  mDataType = mInterface->dataType( mBand );
97 
98  switch ( mDataType )
99  {
105  mIsInteger = true;
106  if ( mDecimalPrecision > -1 )
107  feedback->reportError( QObject::tr( "Input raster is of byte or integer type. The cell values cannot be rounded and will be output using the same data type." ), false );
108  break;
109  default:
110  mIsInteger = false;
111  break;
112  }
113 
114  mInputNoDataValue = inputRaster->dataProvider()->sourceNoDataValue( mBand );
115  mExtent = inputRaster->extent();
116  mLayerWidth = inputRaster->width();
117  mLayerHeight = inputRaster->height();
118  mCrs = inputRaster->crs();
119  mNbCellsXProvider = mInterface->xSize();
120  mNbCellsYProvider = mInterface->ySize();
121  return true;
122 }
123 
124 QVariantMap QgsRoundRasterValuesAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
125 {
126  //prepare output dataset
127  const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context );
128  const QFileInfo fi( outputFile );
129  const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
130  std::unique_ptr< QgsRasterFileWriter > writer = std::make_unique< QgsRasterFileWriter >( outputFile );
131  writer->setOutputProviderKey( QStringLiteral( "gdal" ) );
132  writer->setOutputFormat( outputFormat );
133  std::unique_ptr<QgsRasterDataProvider > provider( writer->createOneBandRaster( mInterface->dataType( mBand ), mNbCellsXProvider, mNbCellsYProvider, mExtent, mCrs ) );
134  if ( !provider )
135  throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( outputFile ) );
136  if ( !provider->isValid() )
137  throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( outputFile, provider->error().message( QgsErrorMessage::Text ) ) );
138 
139  //prepare output provider
140  QgsRasterDataProvider *destinationRasterProvider;
141  destinationRasterProvider = provider.get();
142  destinationRasterProvider->setEditable( true );
143  destinationRasterProvider->setNoDataValue( 1, mInputNoDataValue );
144 
147  const int nbBlocksWidth = static_cast< int >( std::ceil( 1.0 * mLayerWidth / maxWidth ) );
148  const int nbBlocksHeight = static_cast< int >( std::ceil( 1.0 * mLayerHeight / maxHeight ) );
149  const int nbBlocks = nbBlocksWidth * nbBlocksHeight;
150 
151  QgsRasterIterator iter( mInterface.get() );
152  iter.startRasterRead( mBand, mLayerWidth, mLayerHeight, mExtent );
153  int iterLeft = 0;
154  int iterTop = 0;
155  int iterCols = 0;
156  int iterRows = 0;
157  std::unique_ptr< QgsRasterBlock > analysisRasterBlock;
158  while ( iter.readNextRasterPart( mBand, iterCols, iterRows, analysisRasterBlock, iterLeft, iterTop ) )
159  {
160  if ( feedback )
161  feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
162  if ( mIsInteger && mDecimalPrecision > -1 )
163  {
164  //nothing to round, just write raster block
165  analysisRasterBlock->setNoDataValue( mInputNoDataValue );
166  destinationRasterProvider->writeBlock( analysisRasterBlock.get(), mBand, iterLeft, iterTop );
167  }
168  else
169  {
170  for ( int row = 0; row < iterRows; row++ )
171  {
172  if ( feedback && feedback->isCanceled() )
173  break;
174  for ( int column = 0; column < iterCols; column++ )
175  {
176  bool isNoData = false;
177  const double val = analysisRasterBlock->valueAndNoData( row, column, isNoData );
178  if ( isNoData )
179  {
180  analysisRasterBlock->setValue( row, column, mInputNoDataValue );
181  }
182  else
183  {
184  double roundedVal = mInputNoDataValue;
185  if ( mRoundingDirection == 0 && mDecimalPrecision < 0 )
186  {
187  roundedVal = roundUpBaseN( val );
188  }
189  else if ( mRoundingDirection == 0 && mDecimalPrecision > -1 )
190  {
191  const double m = ( val < 0.0 ) ? -1.0 : 1.0;
192  roundedVal = roundUp( val, m );
193  }
194  else if ( mRoundingDirection == 1 && mDecimalPrecision < 0 )
195  {
196  roundedVal = roundNearestBaseN( val );
197  }
198  else if ( mRoundingDirection == 1 && mDecimalPrecision > -1 )
199  {
200  const double m = ( val < 0.0 ) ? -1.0 : 1.0;
201  roundedVal = roundNearest( val, m );
202  }
203  else if ( mRoundingDirection == 2 && mDecimalPrecision < 0 )
204  {
205  roundedVal = roundDownBaseN( val );
206  }
207  else
208  {
209  const double m = ( val < 0.0 ) ? -1.0 : 1.0;
210  roundedVal = roundDown( val, m );
211  }
212  //integer values get automatically cast to double when reading and back to int when writing
213  analysisRasterBlock->setValue( row, column, roundedVal );
214  }
215  }
216  }
217  destinationRasterProvider->writeBlock( analysisRasterBlock.get(), mBand, iterLeft, iterTop );
218  }
219  }
220  destinationRasterProvider->setEditable( false );
221 
222  QVariantMap outputs;
223  outputs.insert( QStringLiteral( "OUTPUT" ), outputFile );
224  return outputs;
225 }
226 
227 double QgsRoundRasterValuesAlgorithm::roundNearest( double value, double m )
228 {
229  return ( std::round( value * m * mScaleFactor ) / mScaleFactor ) * m;
230 }
231 
232 double QgsRoundRasterValuesAlgorithm::roundUp( double value, double m )
233 {
234  return ( std::ceil( value * m * mScaleFactor ) / mScaleFactor ) * m;
235 }
236 
237 double QgsRoundRasterValuesAlgorithm::roundDown( double value, double m )
238 {
239  return ( std::floor( value * m * mScaleFactor ) / mScaleFactor ) * m;
240 }
241 
242 
243 double QgsRoundRasterValuesAlgorithm::roundNearestBaseN( double value )
244 {
245  return static_cast<double>( mMultipleOfBaseN * round( value / mMultipleOfBaseN ) );
246 }
247 
248 double QgsRoundRasterValuesAlgorithm::roundUpBaseN( double value )
249 {
250  return static_cast<double>( mMultipleOfBaseN * ceil( value / mMultipleOfBaseN ) );
251 }
252 
253 double QgsRoundRasterValuesAlgorithm::roundDownBaseN( double value )
254 {
255  return static_cast<double>( mMultipleOfBaseN * floor( value / mMultipleOfBaseN ) );
256 }
257 
258 
259 
260 
@ Int16
Sixteen bit signed integer (qint16)
@ UInt16
Sixteen bit unsigned integer (quint16)
@ Byte
Eight bit unsigned integer (quint8)
@ Int32
Thirty two bit signed integer (qint32)
@ UInt32
Thirty two bit unsigned integer (quint32)
bool isCanceled() const SIP_HOLDGIL
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
virtual QgsRectangle extent() const
Returns the extent of the layer.
QgsCoordinateReferenceSystem crs
Definition: qgsmaplayer.h:79
Contains information about the context in which a processing algorithm is executed.
Custom exception class for processing related exceptions.
Definition: qgsexception.h:83
Base class for providing feedback from a processing algorithm.
virtual void reportError(const QString &error, bool fatalError=false)
Reports that the algorithm encountered an error while executing.
A raster band parameter for Processing algorithms.
@ FlagAdvanced
Parameter is an advanced parameter which should be hidden from users by default.
An enum based parameter for processing algorithms, allowing for selection from predefined values.
A numeric parameter for processing algorithms.
A raster layer destination parameter, for specifying the destination path for a raster layer created ...
A raster layer parameter for processing algorithms.
Base class for raster data providers.
virtual bool setNoDataValue(int bandNo, double noDataValue)
Set no data value on created dataset.
virtual double sourceNoDataValue(int bandNo) const
Value representing no data value.
QgsRasterDataProvider * clone() const override=0
Clone itself, create deep copy.
bool writeBlock(QgsRasterBlock *block, int band, int xOffset=0, int yOffset=0)
Writes pixel data from a raster block into the provider data source.
virtual bool setEditable(bool enabled)
Turns on/off editing mode of the provider.
static QString driverForExtension(const QString &extension)
Returns the GDAL driver name for a specified file extension.
Iterator for sequentially processing raster cells.
static const int DEFAULT_MAXIMUM_TILE_WIDTH
Default maximum tile width.
static const int DEFAULT_MAXIMUM_TILE_HEIGHT
Default maximum tile height.
Represents a raster layer.
int height() const
Returns the height of the (unclipped) raster.
int bandCount() const
Returns the number of bands in this layer.
QgsRasterDataProvider * dataProvider() override
Returns the source data provider.
int width() const
Returns the width of the (unclipped) raster.