Skip to content
Permalink
Browse files

[processing] rescale raster algorithm (fix #26099)

Rescales raster values to the new min-max range preserving shape of the
raster's histogram.
  • Loading branch information
alexbruy committed Jul 8, 2020
1 parent 1878fc7 commit f8d3ac7e3a39bbec69ad8cff9f0b1451c8175357
@@ -142,6 +142,7 @@ SET(QGIS_ANALYSIS_SRCS
processing/qgsalgorithmrenamelayer.cpp
processing/qgsalgorithmrenametablefield.cpp
processing/qgsalgorithmrepairshapefile.cpp
processing/qgsalgorithmrescaleraster.cpp
processing/qgsalgorithmreverselinedirection.cpp
processing/qgsalgorithmrotate.cpp
processing/qgsalgorithmroundrastervalues.cpp
@@ -0,0 +1,168 @@
/***************************************************************************
qgsalgorithmrescaleraster.cpp
---------------------
begin : July 2020
copyright : (C) 2020 by Alexander Bruy
email : alexander dot bruy at gmail dot com
***************************************************************************/

/***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************/

#include <limits>
#include "math.h"
#include "qgsalgorithmrescaleraster.h"
#include "qgsrasterfilewriter.h"

///@cond PRIVATE

QString QgsRescaleRasterAlgorithm::name() const
{
return QStringLiteral( "rescaleraster" );
}

QString QgsRescaleRasterAlgorithm::displayName() const
{
return QObject::tr( "Rescale raster" );
}

QStringList QgsRescaleRasterAlgorithm::tags() const
{
return QObject::tr( "raster,rescale,min,max" ).split( ',' );
}

QString QgsRescaleRasterAlgorithm::group() const
{
return QObject::tr( "Raster analysis" );
}

QString QgsRescaleRasterAlgorithm::groupId() const
{
return QStringLiteral( "rasteranalysis" );
}

QString QgsRescaleRasterAlgorithm::shortHelpString() const
{
return QObject::tr( "Rescales raster layer to the new values range preserving shape "
"(distribution) of the raster's histogram (pixel values)." );
}

QgsRescaleRasterAlgorithm *QgsRescaleRasterAlgorithm::createInstance() const
{
return new QgsRescaleRasterAlgorithm();
}

void QgsRescaleRasterAlgorithm::initAlgorithm( const QVariantMap & )
{
addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT" ), QStringLiteral( "Input raster" ) ) );
addParameter( new QgsProcessingParameterBand( QStringLiteral( "BAND" ), QObject::tr( "Band number" ), 1, QStringLiteral( "INPUT" ) ) );
addParameter( new QgsProcessingParameterNumber( QStringLiteral( "MINIMUM" ), QObject::tr( "New minimum value" ), QgsProcessingParameterNumber::Double, 0 ) );
addParameter( new QgsProcessingParameterNumber( QStringLiteral( "MAXIMUM" ), QObject::tr( "New maximum value" ), QgsProcessingParameterNumber::Double, 255 ) );
addParameter( new QgsProcessingParameterRasterDestination( QStringLiteral( "OUTPUT" ), QObject::tr( "Rescaled" ) ) );
}

bool QgsRescaleRasterAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
{
Q_UNUSED( feedback );

QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT" ), context );
if ( !layer )
throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT" ) ) );

mBand = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
if ( mBand < 1 || mBand > layer->bandCount() )
throw QgsProcessingException( QObject::tr( "Invalid band number for BAND (%1): Valid values for input raster are 1 to %2" )
.arg( mBand )
.arg( layer->bandCount() ) );

mMinimum = parameterAsDouble( parameters, QStringLiteral( "MINIMUM" ), context );
mMaximum = parameterAsDouble( parameters, QStringLiteral( "MAXIMUM" ), context );

mInterface.reset( layer->dataProvider()->clone() );

mCrs = layer->crs();
mLayerWidth = layer->width();
mLayerHeight = layer->height();
mExtent = layer->extent();
mNoData = layer->dataProvider()->sourceNoDataValue( mBand );
mXSize = mInterface->xSize();
mYSize = mInterface->ySize();

return true;
}

QVariantMap QgsRescaleRasterAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
{
feedback->pushInfo( QObject::tr( "Calculating raster minimum and maximum values…" ) );
QgsRasterBandStats stats = mInterface->bandStatistics( mBand, QgsRasterBandStats::Min | QgsRasterBandStats::Max, QgsRectangle(), 0 );

const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context );
QFileInfo fi( outputFile );
const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
std::unique_ptr< QgsRasterFileWriter > writer = qgis::make_unique< QgsRasterFileWriter >( outputFile );
writer->setOutputProviderKey( QStringLiteral( "gdal" ) );
writer->setOutputFormat( outputFormat );
std::unique_ptr<QgsRasterDataProvider > provider( writer->createOneBandRaster( Qgis::Float32, mXSize, mYSize, mExtent, mCrs ) );
if ( !provider )
throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( outputFile ) );
if ( !provider->isValid() )
throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( outputFile, provider->error().message( QgsErrorMessage::Text ) ) );

QgsRasterDataProvider *destProvider = provider.get();
destProvider->setEditable( true );
destProvider->setNoDataValue( 1, mNoData );

int blockWidth = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_WIDTH;
int blockHeight = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_HEIGHT;
int numBlocksX = static_cast< int >( std::ceil( 1.0 * mLayerWidth / blockWidth ) );
int numBlocksY = static_cast< int >( std::ceil( 1.0 * mLayerHeight / blockHeight ) );
int numBlocks = numBlocksX * numBlocksY;

QgsRasterIterator iter( mInterface.get() );
iter.startRasterRead( mBand, mLayerWidth, mLayerHeight, mExtent );
int iterLeft = 0;
int iterTop = 0;
int iterCols = 0;
int iterRows = 0;
std::unique_ptr< QgsRasterBlock > block;
while ( iter.readNextRasterPart( mBand, iterCols, iterRows, block, iterLeft, iterTop ) )
{
if ( feedback )
feedback->setProgress( 100 * ( ( iterTop / blockHeight * numBlocksX ) + iterLeft / blockWidth ) / numBlocks );

for ( int row = 0; row < iterRows; row++ )
{
if ( feedback && feedback->isCanceled() )
break;

for ( int col = 0; col < iterCols; col++ )
{
bool isNoData = false;
double val = block->valueAndNoData( row, col, isNoData );
if ( isNoData )
{
block->setValue( row, col, mNoData );
}
else
{
double newValue = ( ( val - stats.minimumValue ) * ( mMaximum - mMinimum ) / ( stats.maximumValue - stats.minimumValue ) ) + mMinimum;
block->setValue( row, col, newValue );
}
}
}
destProvider->writeBlock( block.get(), mBand, iterLeft, iterTop );
}
destProvider->setEditable( false );

QVariantMap outputs;
outputs.insert( QStringLiteral( "OUTPUT" ), outputFile );
return outputs;
}

///@endcond
@@ -0,0 +1,72 @@
/***************************************************************************
qgsalgorithmrescaleraster.h
------------------------------
begin : July 2020
copyright : (C) 2020 by Alexander Bruy
email : alexander dot bruy at gmail dot com
***************************************************************************/

/***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************/

#ifndef QGSALGORITHMRESCALERASTER_H
#define QGSALGORITHMRESCALERASTER_H

#define SIP_NO_FILE

#include "qgis_sip.h"
#include "qgsprocessingalgorithm.h"
#include "qgsapplication.h"

///@cond PRIVATE

/**
* Native rescale raster algorithm.
*/
class QgsRescaleRasterAlgorithm : public QgsProcessingAlgorithm
{

public:

QgsRescaleRasterAlgorithm() = default;
void initAlgorithm( const QVariantMap &configuration = QVariantMap() ) override;
QString name() const override;
QString displayName() const override;
QStringList tags() const override;
QString group() const override;
QString groupId() const override;
QString shortHelpString() const override;
QgsRescaleRasterAlgorithm *createInstance() const override SIP_FACTORY;

protected:

bool prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;
QVariantMap processAlgorithm( const QVariantMap &parameters,
QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;

private:

int mBand = 1;
int mLayerWidth = 0;
int mLayerHeight = 0;
int mXSize = 0;
int mYSize = 0;
double mNoData;
double mMinimum;
double mMaximum;

Qgis::DataType mDataType;
QgsRectangle mExtent;
QgsCoordinateReferenceSystem mCrs;
std::unique_ptr< QgsRasterInterface > mInterface;
};

///@endcond PRIVATE

#endif // QGSALGORITHMRESCALERASTER_H
@@ -137,6 +137,7 @@
#include "qgsalgorithmrenamelayer.h"
#include "qgsalgorithmrenametablefield.h"
#include "qgsalgorithmrepairshapefile.h"
#include "qgsalgorithmrescaleraster.h"
#include "qgsalgorithmreverselinedirection.h"
#include "qgsalgorithmrotate.h"
#include "qgsalgorithmroundrastervalues.h"
@@ -371,6 +372,7 @@ void QgsNativeAlgorithms::loadAlgorithms()
addAlgorithm( new QgsRenameLayerAlgorithm() );
addAlgorithm( new QgsRenameTableFieldAlgorithm() );
addAlgorithm( new QgsRepairShapefileAlgorithm() );
addAlgorithm( new QgsRescaleRasterAlgorithm() );
addAlgorithm( new QgsReverseLineDirectionAlgorithm() );
addAlgorithm( new QgsRotateFeaturesAlgorithm() );
addAlgorithm( new QgsRoundRasterValuesAlgorithm() );

0 comments on commit f8d3ac7

Please sign in to comment.
You can’t perform that action at this time.