QGIS API Documentation 3.40.0-Bratislava (b56115d8743)
Loading...
Searching...
No Matches
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
25QString QgsRescaleRasterAlgorithm::name() const
26{
27 return QStringLiteral( "rescaleraster" );
28}
29
30QString QgsRescaleRasterAlgorithm::displayName() const
31{
32 return QObject::tr( "Rescale raster" );
33}
34
35QStringList QgsRescaleRasterAlgorithm::tags() const
36{
37 return QObject::tr( "raster,rescale,minimum,maximum,range" ).split( ',' );
38}
39
40QString QgsRescaleRasterAlgorithm::group() const
41{
42 return QObject::tr( "Raster analysis" );
43}
44
45QString QgsRescaleRasterAlgorithm::groupId() const
46{
47 return QStringLiteral( "rasteranalysis" );
48}
49
50QString 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 the original NoData value, but there is "
57 "an option to override it." );
58}
59
60QgsRescaleRasterAlgorithm *QgsRescaleRasterAlgorithm::createInstance() const
61{
62 return new QgsRescaleRasterAlgorithm();
63}
64
65void 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" ), Qgis::ProcessingNumberParameterType::Double, 0 ) );
70 addParameter( new QgsProcessingParameterNumber( QStringLiteral( "MAXIMUM" ), QObject::tr( "New maximum value" ), Qgis::ProcessingNumberParameterType::Double, 255 ) );
71 addParameter( new QgsProcessingParameterNumber( QStringLiteral( "NODATA" ), QObject::tr( "New NoData value" ), Qgis::ProcessingNumberParameterType::Double, QVariant(), true ) );
72
73 std::unique_ptr< QgsProcessingParameterString > createOptsParam = std::make_unique< QgsProcessingParameterString >( QStringLiteral( "CREATE_OPTIONS" ), QObject::tr( "Creation options" ), QVariant(), false, true );
74 createOptsParam->setMetadata( QVariantMap( {{QStringLiteral( "widget_wrapper" ), QVariantMap( {{QStringLiteral( "widget_type" ), QStringLiteral( "rasteroptions" ) }} ) }} ) );
75 createOptsParam->setFlags( createOptsParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
76 addParameter( createOptsParam.release() );
77
78 addParameter( new QgsProcessingParameterRasterDestination( QStringLiteral( "OUTPUT" ), QObject::tr( "Rescaled" ) ) );
79}
80
81bool QgsRescaleRasterAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
82{
83 Q_UNUSED( feedback );
84
85 QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT" ), context );
86 if ( !layer )
87 throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT" ) ) );
88
89 mBand = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
90 if ( mBand < 1 || mBand > layer->bandCount() )
91 throw QgsProcessingException( QObject::tr( "Invalid band number for BAND (%1): Valid values for input raster are 1 to %2" )
92 .arg( mBand )
93 .arg( layer->bandCount() ) );
94
95 mMinimum = parameterAsDouble( parameters, QStringLiteral( "MINIMUM" ), context );
96 mMaximum = parameterAsDouble( parameters, QStringLiteral( "MAXIMUM" ), context );
97
98 mInterface.reset( layer->dataProvider()->clone() );
99
100 mCrs = layer->crs();
101 mLayerWidth = layer->width();
102 mLayerHeight = layer->height();
103 mExtent = layer->extent();
104 if ( parameters.value( QStringLiteral( "NODATA" ) ).isValid() )
105 {
106 mNoData = parameterAsDouble( parameters, QStringLiteral( "NODATA" ), context );
107 }
108 else
109 {
110 mNoData = layer->dataProvider()->sourceNoDataValue( mBand );
111 }
112
113 if ( std::isfinite( mNoData ) )
114 {
115 // Clamp nodata to float32 range, since that's the type of the raster
116 if ( mNoData < std::numeric_limits<float>::lowest() )
117 mNoData = std::numeric_limits<float>::lowest();
118 else if ( mNoData > std::numeric_limits<float>::max() )
119 mNoData = std::numeric_limits<float>::max();
120 }
121
122 mXSize = mInterface->xSize();
123 mYSize = mInterface->ySize();
124
125 return true;
126}
127
128QVariantMap QgsRescaleRasterAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
129{
130 feedback->pushInfo( QObject::tr( "Calculating raster minimum and maximum values…" ) );
131 const QgsRasterBandStats stats = mInterface->bandStatistics( mBand, Qgis::RasterBandStatistic::Min | Qgis::RasterBandStatistic::Max, QgsRectangle(), 0 );
132
133 feedback->pushInfo( QObject::tr( "Rescaling values…" ) );
134 const QString createOptions = parameterAsString( parameters, QStringLiteral( "CREATE_OPTIONS" ), context ).trimmed();
135 const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context );
136 const QFileInfo fi( outputFile );
137 const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
138 std::unique_ptr< QgsRasterFileWriter > writer = std::make_unique< QgsRasterFileWriter >( outputFile );
139 writer->setOutputProviderKey( QStringLiteral( "gdal" ) );
140 if ( !createOptions.isEmpty() )
141 {
142 writer->setCreateOptions( createOptions.split( '|' ) );
143 }
144
145 writer->setOutputFormat( outputFormat );
146 std::unique_ptr<QgsRasterDataProvider > provider( writer->createOneBandRaster( Qgis::DataType::Float32, mXSize, mYSize, mExtent, mCrs ) );
147 if ( !provider )
148 throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( outputFile ) );
149 if ( !provider->isValid() )
150 throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( outputFile, provider->error().message( QgsErrorMessage::Text ) ) );
151
152 QgsRasterDataProvider *destProvider = provider.get();
153 destProvider->setEditable( true );
154 destProvider->setNoDataValue( 1, mNoData );
155
158 const int numBlocksX = static_cast< int >( std::ceil( 1.0 * mLayerWidth / blockWidth ) );
159 const int numBlocksY = static_cast< int >( std::ceil( 1.0 * mLayerHeight / blockHeight ) );
160 const int numBlocks = numBlocksX * numBlocksY;
161
162 QgsRasterIterator iter( mInterface.get() );
163 iter.startRasterRead( mBand, mLayerWidth, mLayerHeight, mExtent );
164 int iterLeft = 0;
165 int iterTop = 0;
166 int iterCols = 0;
167 int iterRows = 0;
168 std::unique_ptr< QgsRasterBlock > inputBlock;
169 while ( iter.readNextRasterPart( mBand, iterCols, iterRows, inputBlock, iterLeft, iterTop ) )
170 {
171 std::unique_ptr< QgsRasterBlock > outputBlock( new QgsRasterBlock( destProvider->dataType( 1 ), iterCols, iterRows ) );
172 feedback->setProgress( 100 * ( ( iterTop / blockHeight * numBlocksX ) + iterLeft / blockWidth ) / numBlocks );
173
174 for ( int row = 0; row < iterRows; row++ )
175 {
176 if ( feedback->isCanceled() )
177 break;
178
179 for ( int col = 0; col < iterCols; col++ )
180 {
181 bool isNoData = false;
182 const double val = inputBlock->valueAndNoData( row, col, isNoData );
183 if ( isNoData )
184 {
185 outputBlock->setValue( row, col, mNoData );
186 }
187 else
188 {
189 const double newValue = ( ( val - stats.minimumValue ) * ( mMaximum - mMinimum ) / ( stats.maximumValue - stats.minimumValue ) ) + mMinimum;
190 outputBlock->setValue( row, col, newValue );
191 }
192 }
193 }
194 destProvider->writeBlock( outputBlock.get(), mBand, iterLeft, iterTop );
195 }
196 destProvider->setEditable( false );
197
198 QVariantMap outputs;
199 outputs.insert( QStringLiteral( "OUTPUT" ), outputFile );
200 return outputs;
201}
202
@ Float32
Thirty two bit floating point (float)
@ Advanced
Parameter is an advanced parameter which should be hidden from users by default.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition qgsfeedback.h:53
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition qgsfeedback.h:61
virtual QgsRectangle extent() const
Returns the extent of the layer.
QgsCoordinateReferenceSystem crs
Definition qgsmaplayer.h:83
Contains information about the context in which a processing algorithm is executed.
Custom exception class for processing related exceptions.
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.
QgsRasterDataProvider * clone() const override=0
Clone itself, create deep copy.
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.
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.