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"
20 #include "qgsalgorithmrescaleraster.h"
21 #include "qgsrasterfilewriter.h"
22
23 ///@cond PRIVATE
24
name() const25 QString QgsRescaleRasterAlgorithm::name() const
26 {
27 return QStringLiteral( "rescaleraster" );
28 }
29
displayName() const30 QString QgsRescaleRasterAlgorithm::displayName() const
31 {
32 return QObject::tr( "Rescale raster" );
33 }
34
tags() const35 QStringList QgsRescaleRasterAlgorithm::tags() const
36 {
37 return QObject::tr( "raster,rescale,minimum,maximum,range" ).split( ',' );
38 }
39
group() const40 QString QgsRescaleRasterAlgorithm::group() const
41 {
42 return QObject::tr( "Raster analysis" );
43 }
44
groupId() const45 QString QgsRescaleRasterAlgorithm::groupId() const
46 {
47 return QStringLiteral( "rasteranalysis" );
48 }
49
shortHelpString() const50 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
createInstance() const60 QgsRescaleRasterAlgorithm *QgsRescaleRasterAlgorithm::createInstance() const
61 {
62 return new QgsRescaleRasterAlgorithm();
63 }
64
initAlgorithm(const QVariantMap &)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
prepareAlgorithm(const QVariantMap & parameters,QgsProcessingContext & context,QgsProcessingFeedback * feedback)75 bool QgsRescaleRasterAlgorithm::prepareAlgorithm( const QVariantMap ¶meters, 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
processAlgorithm(const QVariantMap & parameters,QgsProcessingContext & context,QgsProcessingFeedback * feedback)122 QVariantMap QgsRescaleRasterAlgorithm::processAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
123 {
124 feedback->pushInfo( QObject::tr( "Calculating raster minimum and maximum values…" ) );
125 const 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 const 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
144 const int blockWidth = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_WIDTH;
145 const int blockHeight = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_HEIGHT;
146 const int numBlocksX = static_cast< int >( std::ceil( 1.0 * mLayerWidth / blockWidth ) );
147 const int numBlocksY = static_cast< int >( std::ceil( 1.0 * mLayerHeight / blockHeight ) );
148 const 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 const double val = inputBlock->valueAndNoData( row, col, isNoData );
171 if ( isNoData )
172 {
173 outputBlock->setValue( row, col, mNoData );
174 }
175 else
176 {
177 const 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
191 ///@endcond
192