Merge pull request #37671 from alexbruy/rescale-raster

Rescale raster algorithm for Processing
qgis · Jul 9, 2020 · 12ddfed · 12ddfed
2 parents 588aff0 + 23859b1
commit 12ddfed
Show file tree

Hide file tree

Showing 5 changed files with 283 additions and 0 deletions.
diff --git a/python/plugins/processing/tests/testdata/qgis_algorithm_tests4.yaml b/python/plugins/processing/tests/testdata/qgis_algorithm_tests4.yaml
@@ -2472,5 +2472,34 @@ tests:
         name: expected/filter_geom_lines_poly.gml
         type: vector
 
+  - algorithm: native:rescaleraster
+    name: Rescale raster (original nodata)
+    params:
+      BAND: 1
+      INPUT:
+        name: dem.tif
+        type: raster
+      MAXIMUM: 65535.0
+      MINIMUM: 0.0
+    results:
+      OUTPUT:
+        hash: 63610fb1f7ec3d6530591b1eb479355b52cedfebf6331cec77f2b35c
+        type: rasterhash
+
+  - algorithm: native:rescaleraster
+    name: Rescale raster (new nodata)
+    params:
+      BAND: 1
+      INPUT:
+        name: dem.tif
+        type: raster
+      MAXIMUM: 65535.0
+      MINIMUM: 0.0
+      NODATA: -255.0
+    results:
+      OUTPUT:
+        hash: 63610fb1f7ec3d6530591b1eb479355b52cedfebf6331cec77f2b35c
+        type: rasterhash
+
 
 # See ../README.md for a description of the file format
diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt
@@ -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

diff --git a/src/analysis/processing/qgsalgorithmrescaleraster.cpp b/src/analysis/processing/qgsalgorithmrescaleraster.cpp
@@ -0,0 +1,180 @@
+/***************************************************************************
+                         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,minimum,maximum,range" ).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 a new value range, while preserving the shape "
+                      "(distribution) of the raster's histogram (pixel values). Input values "
+                      "are mapped using a linear interpolation from the source raster's minimum "
+                      "and maximum pixel values to the destination minimum and maximum pixel range.\n\n"
+                      "By default the algorithm preserves original the NODATA value, but there is "
+                      "an option to override it." );
+}
+
+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 QgsProcessingParameterNumber( QStringLiteral( "NODATA" ), QObject::tr( "New NODATA value" ), QgsProcessingParameterNumber::Double, QVariant(), true ) );
+  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();
+  if ( parameters.value( QStringLiteral( "NODATA" ) ).isValid() )
+  {
+    mNoData = parameterAsDouble( parameters, QStringLiteral( "NODATA" ), context );
+  }
+  else
+  {
+    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 );
+
+  feedback->pushInfo( QObject::tr( "Rescaling values…" ) );
+  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 ) )
+  {
+    feedback->setProgress( 100 * ( ( iterTop / blockHeight * numBlocksX ) + iterLeft / blockWidth ) / numBlocks );
+
+    for ( int row = 0; row < iterRows; row++ )
+    {
+      if ( 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
diff --git a/src/analysis/processing/qgsalgorithmrescaleraster.h b/src/analysis/processing/qgsalgorithmrescaleraster.h
@@ -0,0 +1,71 @@
+/***************************************************************************
+                         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;
+
+    QgsRectangle mExtent;
+    QgsCoordinateReferenceSystem mCrs;
+    std::unique_ptr< QgsRasterInterface > mInterface;
+};
+
+///@endcond PRIVATE
+
+#endif // QGSALGORITHMRESCALERASTER_H
diff --git a/src/analysis/processing/qgsnativealgorithms.cpp b/src/analysis/processing/qgsnativealgorithms.cpp
@@ -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() );