QGIS API Documentation  3.20.0-Odense (decaadbb31)
qgsalgorithmrescaleraster.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsalgorithmrescaleraster.cpp
3  ---------------------
4  begin : July 2020
5  copyright : (C) 2020 by Alexander Bruy
6  email : alexander dot bruy 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 <limits>
19 #include "math.h"
21 #include "qgsrasterfilewriter.h"
22 
24 
25 QString QgsRescaleRasterAlgorithm::name() const
26 {
27  return QStringLiteral( "rescaleraster" );
28 }
29 
30 QString QgsRescaleRasterAlgorithm::displayName() const
31 {
32  return QObject::tr( "Rescale raster" );
33 }
34 
35 QStringList QgsRescaleRasterAlgorithm::tags() const
36 {
37  return QObject::tr( "raster,rescale,minimum,maximum,range" ).split( ',' );
38 }
39 
40 QString QgsRescaleRasterAlgorithm::group() const
41 {
42  return QObject::tr( "Raster analysis" );
43 }
44 
45 QString QgsRescaleRasterAlgorithm::groupId() const
46 {
47  return QStringLiteral( "rasteranalysis" );
48 }
49 
50 QString QgsRescaleRasterAlgorithm::shortHelpString() const
51 {
52  return QObject::tr( "Rescales raster layer to a new value range, while preserving the shape "
53  "(distribution) of the raster's histogram (pixel values). Input values "
54  "are mapped using a linear interpolation from the source raster's minimum "
55  "and maximum pixel values to the destination minimum and maximum pixel range.\n\n"
56  "By default the algorithm preserves original the NODATA value, but there is "
57  "an option to override it." );
58 }
59 
60 QgsRescaleRasterAlgorithm *QgsRescaleRasterAlgorithm::createInstance() const
61 {
62  return new QgsRescaleRasterAlgorithm();
63 }
64 
65 void QgsRescaleRasterAlgorithm::initAlgorithm( const QVariantMap & )
66 {
67  addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT" ), QStringLiteral( "Input raster" ) ) );
68  addParameter( new QgsProcessingParameterBand( QStringLiteral( "BAND" ), QObject::tr( "Band number" ), 1, QStringLiteral( "INPUT" ) ) );
69  addParameter( new QgsProcessingParameterNumber( QStringLiteral( "MINIMUM" ), QObject::tr( "New minimum value" ), QgsProcessingParameterNumber::Double, 0 ) );
70  addParameter( new QgsProcessingParameterNumber( QStringLiteral( "MAXIMUM" ), QObject::tr( "New maximum value" ), QgsProcessingParameterNumber::Double, 255 ) );
71  addParameter( new QgsProcessingParameterNumber( QStringLiteral( "NODATA" ), QObject::tr( "New NODATA value" ), QgsProcessingParameterNumber::Double, QVariant(), true ) );
72  addParameter( new QgsProcessingParameterRasterDestination( QStringLiteral( "OUTPUT" ), QObject::tr( "Rescaled" ) ) );
73 }
74 
75 bool QgsRescaleRasterAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
76 {
77  Q_UNUSED( feedback );
78 
79  QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT" ), context );
80  if ( !layer )
81  throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT" ) ) );
82 
83  mBand = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
84  if ( mBand < 1 || mBand > layer->bandCount() )
85  throw QgsProcessingException( QObject::tr( "Invalid band number for BAND (%1): Valid values for input raster are 1 to %2" )
86  .arg( mBand )
87  .arg( layer->bandCount() ) );
88 
89  mMinimum = parameterAsDouble( parameters, QStringLiteral( "MINIMUM" ), context );
90  mMaximum = parameterAsDouble( parameters, QStringLiteral( "MAXIMUM" ), context );
91 
92  mInterface.reset( layer->dataProvider()->clone() );
93 
94  mCrs = layer->crs();
95  mLayerWidth = layer->width();
96  mLayerHeight = layer->height();
97  mExtent = layer->extent();
98  if ( parameters.value( QStringLiteral( "NODATA" ) ).isValid() )
99  {
100  mNoData = parameterAsDouble( parameters, QStringLiteral( "NODATA" ), context );
101  }
102  else
103  {
104  mNoData = layer->dataProvider()->sourceNoDataValue( mBand );
105  }
106 
107  if ( std::isfinite( mNoData ) )
108  {
109  // Clamp nodata to float32 range, since that's the type of the raster
110  if ( mNoData < std::numeric_limits<float>::lowest() )
111  mNoData = std::numeric_limits<float>::lowest();
112  else if ( mNoData > std::numeric_limits<float>::max() )
113  mNoData = std::numeric_limits<float>::max();
114  }
115 
116  mXSize = mInterface->xSize();
117  mYSize = mInterface->ySize();
118 
119  return true;
120 }
121 
122 QVariantMap QgsRescaleRasterAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
123 {
124  feedback->pushInfo( QObject::tr( "Calculating raster minimum and maximum values…" ) );
125  QgsRasterBandStats stats = mInterface->bandStatistics( mBand, QgsRasterBandStats::Min | QgsRasterBandStats::Max, QgsRectangle(), 0 );
126 
127  feedback->pushInfo( QObject::tr( "Rescaling values…" ) );
128  const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context );
129  QFileInfo fi( outputFile );
130  const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
131  std::unique_ptr< QgsRasterFileWriter > writer = std::make_unique< QgsRasterFileWriter >( outputFile );
132  writer->setOutputProviderKey( QStringLiteral( "gdal" ) );
133  writer->setOutputFormat( outputFormat );
134  std::unique_ptr<QgsRasterDataProvider > provider( writer->createOneBandRaster( Qgis::DataType::Float32, mXSize, mYSize, mExtent, mCrs ) );
135  if ( !provider )
136  throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( outputFile ) );
137  if ( !provider->isValid() )
138  throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( outputFile, provider->error().message( QgsErrorMessage::Text ) ) );
139 
140  QgsRasterDataProvider *destProvider = provider.get();
141  destProvider->setEditable( true );
142  destProvider->setNoDataValue( 1, mNoData );
143 
146  int numBlocksX = static_cast< int >( std::ceil( 1.0 * mLayerWidth / blockWidth ) );
147  int numBlocksY = static_cast< int >( std::ceil( 1.0 * mLayerHeight / blockHeight ) );
148  int numBlocks = numBlocksX * numBlocksY;
149 
150  QgsRasterIterator iter( mInterface.get() );
151  iter.startRasterRead( mBand, mLayerWidth, mLayerHeight, mExtent );
152  int iterLeft = 0;
153  int iterTop = 0;
154  int iterCols = 0;
155  int iterRows = 0;
156  std::unique_ptr< QgsRasterBlock > inputBlock;
157  while ( iter.readNextRasterPart( mBand, iterCols, iterRows, inputBlock, iterLeft, iterTop ) )
158  {
159  std::unique_ptr< QgsRasterBlock > outputBlock( new QgsRasterBlock( destProvider->dataType( 1 ), iterCols, iterRows ) );
160  feedback->setProgress( 100 * ( ( iterTop / blockHeight * numBlocksX ) + iterLeft / blockWidth ) / numBlocks );
161 
162  for ( int row = 0; row < iterRows; row++ )
163  {
164  if ( feedback->isCanceled() )
165  break;
166 
167  for ( int col = 0; col < iterCols; col++ )
168  {
169  bool isNoData = false;
170  double val = inputBlock->valueAndNoData( row, col, isNoData );
171  if ( isNoData )
172  {
173  outputBlock->setValue( row, col, mNoData );
174  }
175  else
176  {
177  double newValue = ( ( val - stats.minimumValue ) * ( mMaximum - mMinimum ) / ( stats.maximumValue - stats.minimumValue ) ) + mMinimum;
178  outputBlock->setValue( row, col, newValue );
179  }
180  }
181  }
182  destProvider->writeBlock( outputBlock.get(), mBand, iterLeft, iterTop );
183  }
184  destProvider->setEditable( false );
185 
186  QVariantMap outputs;
187  outputs.insert( QStringLiteral( "OUTPUT" ), outputFile );
188  return outputs;
189 }
190 
@ Float32
Thirty two bit floating point (float)
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:76
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 pushInfo(const QString &info)
Pushes a general informational message from the algorithm.
A raster band parameter for Processing algorithms.
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.
The RasterBandStats struct is a container for statistics about a single raster band.
double minimumValue
The minimum cell value in the raster band.
double maximumValue
The maximum cell value in the raster band.
Raster data container.
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.
Qgis::DataType dataType(int bandNo) const override=0
Returns data type for the band specified by number.
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.
A rectangle specified with double values.
Definition: qgsrectangle.h:42