Branch data Line data Source code
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 : :
25 : 0 : QString QgsRescaleRasterAlgorithm::name() const
26 : : {
27 : 0 : return QStringLiteral( "rescaleraster" );
28 : : }
29 : :
30 : 0 : QString QgsRescaleRasterAlgorithm::displayName() const
31 : : {
32 : 0 : return QObject::tr( "Rescale raster" );
33 : : }
34 : :
35 : 0 : QStringList QgsRescaleRasterAlgorithm::tags() const
36 : : {
37 : 0 : return QObject::tr( "raster,rescale,minimum,maximum,range" ).split( ',' );
38 : 0 : }
39 : :
40 : 0 : QString QgsRescaleRasterAlgorithm::group() const
41 : : {
42 : 0 : return QObject::tr( "Raster analysis" );
43 : : }
44 : :
45 : 0 : QString QgsRescaleRasterAlgorithm::groupId() const
46 : : {
47 : 0 : return QStringLiteral( "rasteranalysis" );
48 : : }
49 : :
50 : 0 : QString QgsRescaleRasterAlgorithm::shortHelpString() const
51 : : {
52 : 0 : 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 : 0 : QgsRescaleRasterAlgorithm *QgsRescaleRasterAlgorithm::createInstance() const
61 : : {
62 : 0 : return new QgsRescaleRasterAlgorithm();
63 : 0 : }
64 : :
65 : 0 : void QgsRescaleRasterAlgorithm::initAlgorithm( const QVariantMap & )
66 : : {
67 : 0 : addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT" ), QStringLiteral( "Input raster" ) ) );
68 : 0 : addParameter( new QgsProcessingParameterBand( QStringLiteral( "BAND" ), QObject::tr( "Band number" ), 1, QStringLiteral( "INPUT" ) ) );
69 : 0 : addParameter( new QgsProcessingParameterNumber( QStringLiteral( "MINIMUM" ), QObject::tr( "New minimum value" ), QgsProcessingParameterNumber::Double, 0 ) );
70 : 0 : addParameter( new QgsProcessingParameterNumber( QStringLiteral( "MAXIMUM" ), QObject::tr( "New maximum value" ), QgsProcessingParameterNumber::Double, 255 ) );
71 : 0 : addParameter( new QgsProcessingParameterNumber( QStringLiteral( "NODATA" ), QObject::tr( "New NODATA value" ), QgsProcessingParameterNumber::Double, QVariant(), true ) );
72 : 0 : addParameter( new QgsProcessingParameterRasterDestination( QStringLiteral( "OUTPUT" ), QObject::tr( "Rescaled" ) ) );
73 : 0 : }
74 : :
75 : 0 : bool QgsRescaleRasterAlgorithm::prepareAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
76 : : {
77 : : Q_UNUSED( feedback );
78 : :
79 : 0 : QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT" ), context );
80 : 0 : if ( !layer )
81 : 0 : throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT" ) ) );
82 : :
83 : 0 : mBand = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
84 : 0 : if ( mBand < 1 || mBand > layer->bandCount() )
85 : 0 : throw QgsProcessingException( QObject::tr( "Invalid band number for BAND (%1): Valid values for input raster are 1 to %2" )
86 : 0 : .arg( mBand )
87 : 0 : .arg( layer->bandCount() ) );
88 : :
89 : 0 : mMinimum = parameterAsDouble( parameters, QStringLiteral( "MINIMUM" ), context );
90 : 0 : mMaximum = parameterAsDouble( parameters, QStringLiteral( "MAXIMUM" ), context );
91 : :
92 : 0 : mInterface.reset( layer->dataProvider()->clone() );
93 : :
94 : 0 : mCrs = layer->crs();
95 : 0 : mLayerWidth = layer->width();
96 : 0 : mLayerHeight = layer->height();
97 : 0 : mExtent = layer->extent();
98 : 0 : if ( parameters.value( QStringLiteral( "NODATA" ) ).isValid() )
99 : : {
100 : 0 : mNoData = parameterAsDouble( parameters, QStringLiteral( "NODATA" ), context );
101 : 0 : }
102 : : else
103 : : {
104 : 0 : mNoData = layer->dataProvider()->sourceNoDataValue( mBand );
105 : : }
106 : :
107 : 0 : if ( std::isfinite( mNoData ) )
108 : : {
109 : : // Clamp nodata to float32 range, since that's the type of the raster
110 : 0 : if ( mNoData < std::numeric_limits<float>::lowest() )
111 : 0 : mNoData = std::numeric_limits<float>::lowest();
112 : 0 : else if ( mNoData > std::numeric_limits<float>::max() )
113 : 0 : mNoData = std::numeric_limits<float>::max();
114 : 0 : }
115 : :
116 : 0 : mXSize = mInterface->xSize();
117 : 0 : mYSize = mInterface->ySize();
118 : :
119 : 0 : return true;
120 : 0 : }
121 : :
122 : 0 : QVariantMap QgsRescaleRasterAlgorithm::processAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
123 : : {
124 : 0 : feedback->pushInfo( QObject::tr( "Calculating raster minimum and maximum values…" ) );
125 : 0 : QgsRasterBandStats stats = mInterface->bandStatistics( mBand, QgsRasterBandStats::Min | QgsRasterBandStats::Max, QgsRectangle(), 0 );
126 : :
127 : 0 : feedback->pushInfo( QObject::tr( "Rescaling values…" ) );
128 : 0 : const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context );
129 : 0 : QFileInfo fi( outputFile );
130 : 0 : const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
131 : 0 : std::unique_ptr< QgsRasterFileWriter > writer = std::make_unique< QgsRasterFileWriter >( outputFile );
132 : 0 : writer->setOutputProviderKey( QStringLiteral( "gdal" ) );
133 : 0 : writer->setOutputFormat( outputFormat );
134 : 0 : std::unique_ptr<QgsRasterDataProvider > provider( writer->createOneBandRaster( Qgis::Float32, mXSize, mYSize, mExtent, mCrs ) );
135 : 0 : if ( !provider )
136 : 0 : throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( outputFile ) );
137 : 0 : if ( !provider->isValid() )
138 : 0 : throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( outputFile, provider->error().message( QgsErrorMessage::Text ) ) );
139 : :
140 : 0 : QgsRasterDataProvider *destProvider = provider.get();
141 : 0 : destProvider->setEditable( true );
142 : 0 : destProvider->setNoDataValue( 1, mNoData );
143 : :
144 : 0 : int blockWidth = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_WIDTH;
145 : 0 : int blockHeight = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_HEIGHT;
146 : 0 : int numBlocksX = static_cast< int >( std::ceil( 1.0 * mLayerWidth / blockWidth ) );
147 : 0 : int numBlocksY = static_cast< int >( std::ceil( 1.0 * mLayerHeight / blockHeight ) );
148 : 0 : int numBlocks = numBlocksX * numBlocksY;
149 : :
150 : 0 : QgsRasterIterator iter( mInterface.get() );
151 : 0 : iter.startRasterRead( mBand, mLayerWidth, mLayerHeight, mExtent );
152 : 0 : int iterLeft = 0;
153 : 0 : int iterTop = 0;
154 : 0 : int iterCols = 0;
155 : 0 : int iterRows = 0;
156 : 0 : std::unique_ptr< QgsRasterBlock > inputBlock;
157 : 0 : while ( iter.readNextRasterPart( mBand, iterCols, iterRows, inputBlock, iterLeft, iterTop ) )
158 : : {
159 : 0 : std::unique_ptr< QgsRasterBlock > outputBlock( new QgsRasterBlock( destProvider->dataType( 1 ), iterCols, iterRows ) );
160 : 0 : feedback->setProgress( 100 * ( ( iterTop / blockHeight * numBlocksX ) + iterLeft / blockWidth ) / numBlocks );
161 : :
162 : 0 : for ( int row = 0; row < iterRows; row++ )
163 : : {
164 : 0 : if ( feedback->isCanceled() )
165 : 0 : break;
166 : :
167 : 0 : for ( int col = 0; col < iterCols; col++ )
168 : : {
169 : 0 : bool isNoData = false;
170 : 0 : double val = inputBlock->valueAndNoData( row, col, isNoData );
171 : 0 : if ( isNoData )
172 : : {
173 : 0 : outputBlock->setValue( row, col, mNoData );
174 : 0 : }
175 : : else
176 : : {
177 : 0 : double newValue = ( ( val - stats.minimumValue ) * ( mMaximum - mMinimum ) / ( stats.maximumValue - stats.minimumValue ) ) + mMinimum;
178 : 0 : outputBlock->setValue( row, col, newValue );
179 : : }
180 : 0 : }
181 : 0 : }
182 : 0 : destProvider->writeBlock( outputBlock.get(), mBand, iterLeft, iterTop );
183 : 0 : }
184 : 0 : destProvider->setEditable( false );
185 : :
186 : 0 : QVariantMap outputs;
187 : 0 : outputs.insert( QStringLiteral( "OUTPUT" ), outputFile );
188 : 0 : return outputs;
189 : 0 : }
190 : :
191 : : ///@endcond
|
