Processing raster rank algorithm (#60605)

src/analysis/CMakeLists.txt

@@ -193,6 +193,7 @@ set(QGIS_ANALYSIS_SRCS
   processing/qgsalgorithmrasterlogicalop.cpp
   processing/qgsalgorithmrasterminmax.cpp
   processing/qgsalgorithmrasterize.cpp
+  processing/qgsalgorithmrasterrank.cpp
   processing/qgsalgorithmrastersampling.cpp
   processing/qgsalgorithmrasterstackposition.cpp
   processing/qgsalgorithmrasterstatistics.cpp

src/analysis/processing/qgsalgorithmrasterrank.cpp

@@ -0,0 +1,327 @@
+/***************************************************************************
+                         qgsalgorithmrasterrank.cpp
+                         ---------------------
+    begin                : February 2024
+    copyright            : (C) 2024 by Mathieu Pellerin
+    email                : mathieu at opengis dot ch
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   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 "qgsalgorithmrasterrank.h"
+#include "qgsrasterfilewriter.h"
+#include "qgsrasterprojector.h"
+
+///@cond PRIVATE
+
+QString QgsRasterRankAlgorithm::name() const
+{
+  return QStringLiteral( "rasterrank" );
+}
+
+QString QgsRasterRankAlgorithm::displayName() const
+{
+  return QObject::tr( "Raster rank" );
+}
+
+QStringList QgsRasterRankAlgorithm::tags() const
+{
+  return QObject::tr( "raster,rank" ).split( ',' );
+}
+
+QString QgsRasterRankAlgorithm::group() const
+{
+  return QObject::tr( "Raster analysis" );
+}
+
+QString QgsRasterRankAlgorithm::groupId() const
+{
+  return QStringLiteral( "rasteranalysis" );
+}
+
+QString QgsRasterRankAlgorithm::shortHelpString() const
+{
+  return QObject::tr( "Performs a cell-by-cell analysis in which output values match the rank of a "
+                      "sorted list of overlapping cell values from input layers. The output raster "
+                      "will be multi-band if multiple ranks are provided.\n"
+                      "If multiband rasters are used in the data raster stack, the algorithm will always "
+                      "perform the analysis on the first band of the rasters." );
+}
+
+QgsRasterRankAlgorithm *QgsRasterRankAlgorithm::createInstance() const
+{
+  return new QgsRasterRankAlgorithm();
+}
+
+void QgsRasterRankAlgorithm::initAlgorithm( const QVariantMap & )
+{
+  addParameter( new QgsProcessingParameterMultipleLayers( QStringLiteral( "INPUT_RASTERS" ), QObject::tr( "Input raster layers" ), Qgis::ProcessingSourceType::Raster ) );
+  auto ranksParameter = std::make_unique<QgsProcessingParameterString>( QStringLiteral( "RANKS" ), QObject::tr( "Rank(s) (separate multiple ranks using commas)" ), 1 );
+  ranksParameter->setHelp( QObject::tr( "A rank value must be numerical, with multiple ranks separated by commas. The rank will be used to "
+                                        "generate output values from sorted lists of input layers’ cell values. A rank value of 1 will pick "
+                                        "the first value from a given sorted input layers’ cell values list (i.e. the minimum value). "
+                                        "Negative rank values are supported, and will behave like a negative index. A rank value of -2 will "
+                                        "pick the second to last value in sorted input values lists, while a rank value of -1 will pick the "
+                                        "last value (i.e. the maximum value)." ) );
+  addParameter( ranksParameter.release() );
+  addParameter( new QgsProcessingParameterEnum( QStringLiteral( "NODATA_HANDLING" ), QObject::tr( "NoData value handling" ), QStringList() << QObject::tr( "Exclude NoData from values lists" ) << QObject::tr( "Presence of NoData in a values list results in NoData output cell" ), false, 0 ) );
+
+  auto extentParam = std::make_unique<QgsProcessingParameterExtent>( QStringLiteral( "EXTENT" ), QObject::tr( "Output extent" ), QVariant(), true );
+  extentParam->setHelp( QObject::tr( "Extent of the output layer. If not specified, the extent will be the overall extent of all input layers" ) );
+  extentParam->setFlags( extentParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
+  addParameter( extentParam.release() );
+  auto cellSizeParam = std::make_unique<QgsProcessingParameterNumber>( QStringLiteral( "CELL_SIZE" ), QObject::tr( "Output cell size" ), Qgis::ProcessingNumberParameterType::Double, QVariant(), true, 0.0 );
+  cellSizeParam->setHelp( QObject::tr( "Cell size of the output layer. If not specified, the smallest cell size from the input layers will be used" ) );
+  cellSizeParam->setFlags( cellSizeParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
+  addParameter( cellSizeParam.release() );
+  auto crsParam = std::make_unique<QgsProcessingParameterCrs>( QStringLiteral( "CRS" ), QObject::tr( "Output CRS" ), QVariant(), true );
+  crsParam->setHelp( QObject::tr( "CRS of the output layer. If not specified, the CRS of the first input layer will be used" ) );
+  crsParam->setFlags( crsParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
+  addParameter( crsParam.release() );
+
+  addParameter( new QgsProcessingParameterRasterDestination( QStringLiteral( "OUTPUT" ), QObject::tr( "Ranked" ) ) );
+}
+
+bool QgsRasterRankAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
+{
+  const QStringList rankStrings = parameterAsString( parameters, QStringLiteral( "RANKS" ), context ).split( QLatin1String( "," ) );
+  for ( const QString &rankString : rankStrings )
+  {
+    bool ok = false;
+    const int rank = rankString.toInt( &ok );
+    if ( ok && rank != 0 )
+    {
+      mRanks << rank;
+    }
+    else
+    {
+      throw QgsProcessingException( QObject::tr( "Rank values must be integers (found \"%1\")" ).arg( rankString ) );
+    }
+  }
+
+  if ( mRanks.isEmpty() )
+  {
+    throw QgsProcessingException( QObject::tr( "No valid non-zero rank value(s) provided" ) );
+    return false;
+  }
+
+  const QList<QgsMapLayer *> layers = parameterAsLayerList( parameters, QStringLiteral( "INPUT_RASTERS" ), context );
+  for ( const QgsMapLayer *layer : std::as_const( layers ) )
+  {
+    if ( !qobject_cast<const QgsRasterLayer *>( layer ) || !layer->dataProvider() )
+      continue;
+
+    std::unique_ptr<QgsMapLayer> clonedLayer;
+    clonedLayer.reset( layer->clone() );
+    clonedLayer->moveToThread( nullptr );
+    mLayers.push_back( std::move( clonedLayer ) );
+  }
+
+  if ( mLayers.empty() )
+  {
+    throw QgsProcessingException( QObject::tr( "At least one raster input layer must be specified" ) );
+    return false;
+  }
+
+  return true;
+}
+
+
+QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
+{
+  QList<QgsMapLayer *> layers; // Needed for QgsProcessingUtils::combineLayerExtents
+  for ( auto &layer : mLayers )
+  {
+    layer->moveToThread( QThread::currentThread() );
+    layers << layer.get();
+  }
+
+  QgsCoordinateReferenceSystem outputCrs;
+  if ( parameters.value( QStringLiteral( "CRS" ) ).isValid() )
+  {
+    outputCrs = parameterAsCrs( parameters, QStringLiteral( "CRS" ), context );
+  }
+  else
+  {
+    outputCrs = mLayers[0]->crs();
+  }
+
+
+  QgsRasterLayer *templateRasterLayer = qobject_cast<QgsRasterLayer *>( mLayers[0].get() );
+  const Qgis::DataType outputDataType = templateRasterLayer->dataProvider()->dataType( 1 );
+  double outputNoData = 0.0;
+  if ( templateRasterLayer->dataProvider()->sourceHasNoDataValue( 1 ) )
+  {
+    outputNoData = templateRasterLayer->dataProvider()->sourceNoDataValue( 1 );
+  }
+  else
+  {
+    outputNoData = -FLT_MAX;
+  }
+  const bool outputNoDataOverride = parameterAsInt( parameters, QStringLiteral( "NODATA_HANDLING" ), context ) == 1;
+
+  QgsRectangle outputExtent;
+  if ( parameters.value( QStringLiteral( "EXTENT" ) ).isValid() )
+  {
+    outputExtent = parameterAsExtent( parameters, QStringLiteral( "EXTENT" ), context, outputCrs );
+  }
+  else
+  {
+    outputExtent = QgsProcessingUtils::combineLayerExtents( layers, outputCrs, context );
+  }
+
+  double minCellSizeX = 1e9;
+  double minCellSizeY = 1e9;
+  for ( auto &layer : mLayers )
+  {
+    QgsRasterLayer *rasterLayer = qobject_cast<QgsRasterLayer *>( layer.get() );
+
+    QgsRectangle extent = rasterLayer->extent();
+    if ( rasterLayer->crs() != outputCrs )
+    {
+      QgsCoordinateTransform ct( rasterLayer->crs(), outputCrs, context.transformContext() );
+      extent = ct.transformBoundingBox( extent );
+    }
+
+    minCellSizeX = std::min( minCellSizeX, ( extent.xMaximum() - extent.xMinimum() ) / rasterLayer->width() );
+    minCellSizeY = std::min( minCellSizeY, ( extent.yMaximum() - extent.yMinimum() ) / rasterLayer->height() );
+  }
+
+  double outputCellSizeX = parameterAsDouble( parameters, QStringLiteral( "CELL_SIZE" ), context );
+  double outputCellSizeY = outputCellSizeX;
+  if ( outputCellSizeX == 0 )
+  {
+    outputCellSizeX = minCellSizeX;
+    outputCellSizeY = minCellSizeY;
+  }
+
+  qgssize cols = static_cast<qgssize>( std::round( ( outputExtent.xMaximum() - outputExtent.xMinimum() ) / outputCellSizeX ) );
+  qgssize rows = static_cast<qgssize>( std::round( ( outputExtent.yMaximum() - outputExtent.yMinimum() ) / outputCellSizeY ) );
+
+  const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context );
+  const QFileInfo fi( outputFile );
+  const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
+
+  auto writer = std::make_unique<QgsRasterFileWriter>( outputFile );
+  writer->setOutputFormat( outputFormat );
+  std::unique_ptr<QgsRasterDataProvider> provider( writer->createMultiBandRaster( outputDataType, cols, rows, outputExtent, outputCrs, mRanks.size() ) );
+  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 ) ) );
+  provider->setNoDataValue( 1, outputNoData );
+
+  std::map<QString, std::unique_ptr<QgsRasterInterface>> newProjectorInterfaces;
+  std::map<QString, QgsRasterInterface *> inputInterfaces;
+  std::map<QString, std::unique_ptr<QgsRasterBlock>> inputBlocks;
+  std::vector<std::unique_ptr<QgsRasterBlock>> outputBlocks;
+  outputBlocks.resize( mRanks.size() );
+  for ( auto &layer : mLayers )
+  {
+    QgsRasterLayer *rasterLayer = qobject_cast<QgsRasterLayer *>( layer.get() );
+    if ( rasterLayer->crs() != outputCrs )
+    {
+      QgsRasterProjector *projector = new QgsRasterProjector();
+      projector->setCrs( rasterLayer->crs(), outputCrs, context.transformContext() );
+      projector->setInput( rasterLayer->dataProvider() );
+      projector->setPrecision( QgsRasterProjector::Exact );
+      newProjectorInterfaces[rasterLayer->id()].reset( projector );
+      inputInterfaces[rasterLayer->id()] = projector;
+    }
+    else
+    {
+      inputInterfaces[rasterLayer->id()] = rasterLayer->dataProvider();
+    }
+  }
+
+  QgsRasterIterator rasterIterator( inputInterfaces[mLayers.at( 0 )->id()] );
+  rasterIterator.startRasterRead( 1, cols, rows, outputExtent );
+  int blockCount = static_cast<int>( rasterIterator.blockCount() );
+
+  const double step = blockCount > 0 ? 100.0 / blockCount : 0;
+  std::vector<double> inputValues;
+  inputValues.resize( mLayers.size() );
+  for ( int currentBlock = 0; currentBlock < blockCount; currentBlock++ )
+  {
+    if ( feedback->isCanceled() )
+    {
+      break;
+    }
+    feedback->setProgress( currentBlock * step );
+
+    int iterLeft = 0;
+    int iterTop = 0;
+    int iterCols = 0;
+    int iterRows = 0;
+    QgsRectangle blockExtent;
+    rasterIterator.next( 1, iterCols, iterRows, iterLeft, iterTop, blockExtent );
+
+    for ( const auto &inputInterface : inputInterfaces )
+    {
+      inputBlocks[inputInterface.first].reset( inputInterface.second->block( 1, blockExtent, iterCols, iterRows ) );
+    }
+
+    for ( int i = 0; i < mRanks.size(); i++ )
+    {
+      outputBlocks[i].reset( new QgsRasterBlock( outputDataType, iterCols, iterRows ) );
+      outputBlocks[i]->setNoDataValue( outputNoData );
+    }
+
+    for ( int row = 0; row < iterRows; row++ )
+    {
+      for ( int col = 0; col < iterCols; col++ )
+      {
+        int valuesCount = 0;
+        for ( const auto &inputBlock : inputBlocks )
+        {
+          bool isNoData = false;
+          const double value = inputBlock.second->valueAndNoData( row, col, isNoData );
+          if ( !isNoData )
+          {
+            inputValues[valuesCount] = value;
+            valuesCount++;
+          }
+          else if ( outputNoDataOverride )
+          {
+            valuesCount = 0;
+            break;
+          }
+        }
+        std::sort( inputValues.begin(), inputValues.begin() + valuesCount );
+
+        for ( int i = 0; i < mRanks.size(); i++ )
+        {
+          if ( valuesCount >= std::abs( mRanks[i] ) )
+          {
+            outputBlocks[i]->setValue( row, col, inputValues[mRanks[i] > 0 ? mRanks[i] - 1 : valuesCount + mRanks[i]] );
+          }
+          else
+          {
+            outputBlocks[i]->setValue( row, col, outputNoData );
+          }
+        }
+      }
+    }
+
+    for ( int i = 0; i < mRanks.size(); i++ )
+    {
+      if ( !provider->writeBlock( outputBlocks[i].get(), i + 1, iterLeft, iterTop ) )
+      {
+        throw QgsProcessingException( QObject::tr( "Could not write output raster block: %1" ).arg( provider->error().summary() ) );
+      }
+    }
+  }
+
+  QVariantMap outputs;
+  outputs.insert( QStringLiteral( "OUTPUT" ), outputFile );
+  return outputs;
+}
+
+///@endcond

src/analysis/processing/qgsalgorithmrasterrank.h

@@ -0,0 +1,55 @@
+/***************************************************************************
+                         qgsalgorithmrasterrank.h
+                         ---------------------
+    begin                : February 2024
+    copyright            : (C) 2024 by Mathieu Pellerin
+    email                : mathieu at opengis dot ch
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   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 QGSALGORITHMRASTERRANK_H
+#define QGSALGORITHMRASTERRANK_H
+
+#define SIP_NO_FILE
+
+#include "qgis_sip.h"
+#include "qgsprocessingalgorithm.h"
+#include "qgsapplication.h"
+
+///@cond PRIVATE
+
+/**
+ * Native raster rank algorithm.
+ */
+class QgsRasterRankAlgorithm : public QgsProcessingAlgorithm
+{
+  public:
+    QgsRasterRankAlgorithm() = 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;
+    QgsRasterRankAlgorithm *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;
+
+    QList<int> mRanks;
+    std::vector<std::unique_ptr<QgsMapLayer>> mLayers;
+};
+
+///@endcond PRIVATE
+
+#endif // QGSALGORITHMRASTERRANK_H

src/analysis/processing/qgsnativealgorithms.cpp

@@ -175,6 +175,7 @@
 #include "qgsalgorithmrasterlogicalop.h"
 #include "qgsalgorithmrasterminmax.h"
 #include "qgsalgorithmrasterize.h"
+#include "qgsalgorithmrasterrank.h"
 #include "qgsalgorithmrastersampling.h"
 #include "qgsalgorithmrasterstackposition.h"
 #include "qgsalgorithmrasterstatistics.h"
@@ -499,6 +500,7 @@ void QgsNativeAlgorithms::loadAlgorithms()
   addAlgorithm( new QgsRasterizeAlgorithm() );
   addAlgorithm( new QgsRasterPixelsToPointsAlgorithm() );
   addAlgorithm( new QgsRasterPixelsToPolygonsAlgorithm() );
+  addAlgorithm( new QgsRasterRankAlgorithm() );
   addAlgorithm( new QgsRasterSamplingAlgorithm() );
   addAlgorithm( new QgsRasterStackHighestPositionAlgorithm() );
   addAlgorithm( new QgsRasterStackLowestPositionAlgorithm() );

tests/src/analysis/testqgsprocessingalgspt1.cpp

@@ -105,6 +105,9 @@ class TestQgsProcessingAlgsPt1 : public QgsTest
     void createConstantRaster_data();
     void createConstantRaster();
 
+    void rasterRank_data();
+    void rasterRank();
+
     void densifyGeometries_data();
     void densifyGeometries();
 
@@ -1897,6 +1900,113 @@ void TestQgsProcessingAlgsPt1::createConstantRaster()
   }
 }
 
+void TestQgsProcessingAlgsPt1::rasterRank_data()
+{
+  QTest::addColumn<QString>( "expectedRaster" );
+  QTest::addColumn<QString>( "ranks" );
+  QTest::addColumn<int>( "nodataHandling" );
+
+  /*
+   * Testcase 1
+   */
+  QTest::newRow( "testcase 1" )
+    << QStringLiteral( "/rasterRank_testcase1.tif" )
+    << QString::number( 2 )
+    << 0;
+
+  /*
+   * Testcase 2
+   */
+  QTest::newRow( "testcase 2" )
+    << QStringLiteral( "/rasterRank_testcase2.tif" )
+    << QString::number( -2 )
+    << 0;
+
+  /*
+   * Testcase 3
+   */
+  QTest::newRow( "testcase 3" )
+    << QStringLiteral( "/rasterRank_testcase3.tif" )
+    << QString::number( 2 )
+    << 1;
+
+  /*
+   * Testcase 4
+   */
+  QTest::newRow( "testcase 4" )
+    << QStringLiteral( "/rasterRank_testcase4.tif" )
+    << QStringLiteral( "2,-2" )
+    << 0;
+}
+
+void TestQgsProcessingAlgsPt1::rasterRank()
+{
+  QFETCH( QString, expectedRaster );
+  QFETCH( QString, ranks );
+  QFETCH( int, nodataHandling );
+
+  //prepare input params
+  std::unique_ptr<QgsProcessingAlgorithm> alg( QgsApplication::processingRegistry()->createAlgorithmById( QStringLiteral( "native:rasterrank" ) ) );
+
+  const QString testDataPath( TEST_DATA_DIR ); //defined in CMakeLists.txt
+
+  QVariantMap parameters;
+
+  parameters.insert( QStringLiteral( "INPUT_RASTERS" ), QStringList() << testDataPath + "/raster/rank1.tif" << testDataPath + "/raster/rank2.tif" << testDataPath + "/raster/rank3.tif" << testDataPath + "/raster/rank4.tif" );
+  parameters.insert( QStringLiteral( "RANKS" ), ranks );
+  parameters.insert( QStringLiteral( "NODATA_HANDLING" ), nodataHandling );
+  parameters.insert( QStringLiteral( "OUTPUT" ), QgsProcessing::TEMPORARY_OUTPUT );
+
+  bool ok = false;
+  QgsProcessingFeedback feedback;
+  QVariantMap results;
+
+  //prepare expected raster
+  auto expectedRasterLayer = std::make_unique<QgsRasterLayer>( testDataPath + "/control_images/expected_rasterRank" + expectedRaster, "expectedDataset", "gdal" );
+  std::unique_ptr<QgsRasterInterface> expectedInterface( expectedRasterLayer->dataProvider()->clone() );
+  QgsRasterIterator expectedIter( expectedInterface.get() );
+  expectedIter.startRasterRead( 1, expectedRasterLayer->width(), expectedRasterLayer->height(), expectedInterface->extent() );
+
+  //run algorithm...
+  auto context = std::make_unique<QgsProcessingContext>();
+  results = alg->run( parameters, *context, &feedback, &ok );
+  QVERIFY( ok );
+
+  //...and check results with expected datasets
+  auto outputRaster = std::make_unique<QgsRasterLayer>( results.value( QStringLiteral( "OUTPUT" ) ).toString(), "output", "gdal" );
+  std::unique_ptr<QgsRasterInterface> outputInterface( outputRaster->dataProvider()->clone() );
+
+  QCOMPARE( outputRaster->width(), expectedRasterLayer->width() );
+  QCOMPARE( outputRaster->height(), expectedRasterLayer->height() );
+
+  QgsRasterIterator outputIter( outputInterface.get() );
+  outputIter.startRasterRead( 1, outputRaster->width(), outputRaster->height(), outputInterface->extent() );
+  int outputIterLeft = 0;
+  int outputIterTop = 0;
+  int outputIterCols = 0;
+  int outputIterRows = 0;
+  int expectedIterLeft = 0;
+  int expectedIterTop = 0;
+  int expectedIterCols = 0;
+  int expectedIterRows = 0;
+
+  std::unique_ptr<QgsRasterBlock> outputRasterBlock;
+  std::unique_ptr<QgsRasterBlock> expectedRasterBlock;
+
+  while ( outputIter.readNextRasterPart( 1, outputIterCols, outputIterRows, outputRasterBlock, outputIterLeft, outputIterTop ) && expectedIter.readNextRasterPart( 1, expectedIterCols, expectedIterRows, expectedRasterBlock, expectedIterLeft, expectedIterTop ) )
+  {
+    for ( int row = 0; row < expectedIterRows; row++ )
+    {
+      for ( int column = 0; column < expectedIterCols; column++ )
+      {
+        const double expectedValue = expectedRasterBlock->value( row, column );
+        const double outputValue = outputRasterBlock->value( row, column );
+        QCOMPARE( outputValue, expectedValue );
+      }
+    }
+  }
+}
+
 
 void TestQgsProcessingAlgsPt1::densifyGeometries_data()
 {