From 2880bced9d2afd5e82a100a17175af029efd75fa Mon Sep 17 00:00:00 2001 From: Mathieu Pellerin Date: Thu, 13 Feb 2025 18:52:39 +0700 Subject: [PATCH 1/9] [processing] Brand new raster rank algorithm --- src/analysis/CMakeLists.txt | 1 + .../processing/qgsalgorithmrasterrank.cpp | 246 ++++++++++++++++++ .../processing/qgsalgorithmrasterrank.h | 54 ++++ .../processing/qgsnativealgorithms.cpp | 2 + 4 files changed, 303 insertions(+) create mode 100644 src/analysis/processing/qgsalgorithmrasterrank.cpp create mode 100644 src/analysis/processing/qgsalgorithmrasterrank.h diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt index 751f0ee52cd9..7c825293b3af 100644 --- a/src/analysis/CMakeLists.txt +++ b/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 diff --git a/src/analysis/processing/qgsalgorithmrasterrank.cpp b/src/analysis/processing/qgsalgorithmrasterrank.cpp new file mode 100644 index 000000000000..c180c87b88ca --- /dev/null +++ b/src/analysis/processing/qgsalgorithmrasterrank.cpp @@ -0,0 +1,246 @@ +/*************************************************************************** + 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( "Performing a cell-by-cell analysis in which output values match the rank of a sorted list of overlapping cell values from input layers." ); +} + +QgsRasterRankAlgorithm *QgsRasterRankAlgorithm::createInstance() const +{ + return new QgsRasterRankAlgorithm(); +} + +void QgsRasterRankAlgorithm::initAlgorithm( const QVariantMap & ) +{ + addParameter( new QgsProcessingParameterMultipleLayers( QStringLiteral( "LAYERS" ), QObject::tr( "Input layers" ), Qgis::ProcessingSourceType::Raster ) ); + addParameter( new QgsProcessingParameterNumber( QStringLiteral( "RANK" ), QObject::tr( "Rank" ), Qgis::ProcessingNumberParameterType::Integer, 1 ) ); + + auto extentParam = std::make_unique( 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( 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( 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( "Calculated" ) ) ); +} + +bool QgsRasterRankAlgorithm::prepareAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) +{ + const QList layers = parameterAsLayerList( parameters, QStringLiteral( "LAYERS" ), context ); + + for ( const QgsMapLayer *layer : std::as_const( layers ) ) + { + if ( !qobject_cast( layer ) || !layer->dataProvider() ) + continue; + + QgsMapLayer *clonedLayer = layer->clone(); + clonedLayer->moveToThread( nullptr ); + mLayers << clonedLayer; + } + + if ( mLayers.isEmpty() ) + { + feedback->reportError( QObject::tr( "No raster layers selected" ), false ); + return false; + } + + return true; +} + + +QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) +{ + for ( QgsMapLayer *layer : std::as_const( mLayers ) ) + { + layer->moveToThread( QThread::currentThread() ); + } + + const int rank = parameterAsInt( parameters, QStringLiteral( "RANK" ), context ); + QgsCoordinateReferenceSystem outputCrs; + if ( parameters.value( QStringLiteral( "CRS" ) ).isValid() ) + { + outputCrs = parameterAsCrs( parameters, QStringLiteral( "CRS" ), context ); + } + else + { + outputCrs = mLayers.at( 0 )->crs(); + } + const Qgis::DataType outputDataType = qobject_cast( mLayers.at( 0 ) )->dataProvider()->dataType( 1 ); + QgsRectangle outputExtent; + if ( parameters.value( QStringLiteral( "EXTENT" ) ).isValid() ) + { + outputExtent = parameterAsExtent( parameters, QStringLiteral( "EXTENT" ), context, outputCrs ); + } + else + { + outputExtent = QgsProcessingUtils::combineLayerExtents( mLayers, outputCrs, context ); + } + + double minCellSize = 1e9; + std::map> inputBlocks; + for ( QgsMapLayer *layer : mLayers ) + { + QgsRasterLayer *rLayer = qobject_cast( layer ); + + QgsRectangle extent = rLayer->extent(); + if ( rLayer->crs() != outputCrs ) + { + QgsCoordinateTransform ct( rLayer->crs(), outputCrs, context.transformContext() ); + extent = ct.transformBoundingBox( extent ); + } + + double cellSize = ( extent.xMaximum() - extent.xMinimum() ) / rLayer->width(); + if ( cellSize < minCellSize ) + { + minCellSize = cellSize; + } + + inputBlocks[rLayer->id()] = std::make_unique(); + } + + double outputCellSize = parameterAsDouble( parameters, QStringLiteral( "CELL_SIZE" ), context ); + if ( outputCellSize == 0 ) + { + outputCellSize = minCellSize; + } + + double cols = std::round( ( outputExtent.xMaximum() - outputExtent.xMinimum() ) / outputCellSize ); + double rows = std::round( ( outputExtent.yMaximum() - outputExtent.yMinimum() ) / outputCellSize ); + + const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context ); + const QFileInfo fi( outputFile ); + const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() ); + + auto writer = std::make_unique( outputFile ); + writer->setOutputFormat( outputFormat ); + std::unique_ptr provider( writer->createOneBandRaster( outputDataType, cols, rows, outputExtent, outputCrs ) ); + 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, 0 ); + + const double step = rows > 0 ? 100.0 / rows : 1; + for ( int row = 0; row < rows; row++ ) + { + if ( feedback->isCanceled() ) + { + break; + } + + QgsRasterBlock block( outputDataType, cols, 1 ); + block.setNoDataValue( 0 ); + + // Calculates the rect for a single row read + QgsRectangle rowExtent( outputExtent ); + rowExtent.setYMaximum( rowExtent.yMaximum() - outputCellSize * row ); + rowExtent.setYMinimum( rowExtent.yMaximum() - outputCellSize ); + + for ( QgsMapLayer *layer : mLayers ) + { + QgsRasterLayer *rLayer = qobject_cast( layer ); + + if ( rLayer->crs() != outputCrs ) + { + QgsRasterProjector proj; + proj.setCrs( rLayer->crs(), outputCrs, context.transformContext() ); + proj.setInput( rLayer->dataProvider() ); + proj.setPrecision( QgsRasterProjector::Exact ); + inputBlocks[rLayer->id()].reset( proj.block( 1, rowExtent, cols, 1 ) ); + } + else + { + inputBlocks[rLayer->id()].reset( rLayer->dataProvider()->block( 1, rowExtent, cols, 1 ) ); + } + } + + for ( int col = 0; col < cols; col++ ) + { + QList values; + for ( const auto &inputBlock : inputBlocks ) + { + bool isNoData = false; + const double value = inputBlock.second->valueAndNoData( 0, col, isNoData ); + if ( !isNoData ) + { + values << value; + } + } + + std::sort( values.begin(), values.end() ); + if ( values.size() >= rank ) + { + block.setValue( 0, col, values.at( rank - 1 ) ); + } + else + { + block.setValue( 0, col, 0 ); + } + } + + provider->writeBlock( &block, 1, 0, row ); + feedback->setProgress( row * step ); + } + + qDeleteAll( mLayers ); + mLayers.clear(); + + QVariantMap outputs; + outputs.insert( QStringLiteral( "OUTPUT" ), outputFile ); + return outputs; +} + +///@endcond diff --git a/src/analysis/processing/qgsalgorithmrasterrank.h b/src/analysis/processing/qgsalgorithmrasterrank.h new file mode 100644 index 000000000000..6b298ace9516 --- /dev/null +++ b/src/analysis/processing/qgsalgorithmrasterrank.h @@ -0,0 +1,54 @@ +/*************************************************************************** + 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 ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override; + QVariantMap processAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override; + + QList mLayers; +}; + +///@endcond PRIVATE + +#endif // QGSALGORITHMRASTERRANK_H diff --git a/src/analysis/processing/qgsnativealgorithms.cpp b/src/analysis/processing/qgsnativealgorithms.cpp index c730592528f5..06e58cc3eb64 100644 --- a/src/analysis/processing/qgsnativealgorithms.cpp +++ b/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() ); From 8df0f889929ccfac696e2e135bcd1c860bbab022 Mon Sep 17 00:00:00 2001 From: Mathieu Pellerin Date: Thu, 13 Feb 2025 19:10:30 +0700 Subject: [PATCH 2/9] Proper handling of no data values --- .../processing/qgsalgorithmrasterrank.cpp | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/analysis/processing/qgsalgorithmrasterrank.cpp b/src/analysis/processing/qgsalgorithmrasterrank.cpp index c180c87b88ca..191df6d38b99 100644 --- a/src/analysis/processing/qgsalgorithmrasterrank.cpp +++ b/src/analysis/processing/qgsalgorithmrasterrank.cpp @@ -118,7 +118,19 @@ QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶met { outputCrs = mLayers.at( 0 )->crs(); } + + const Qgis::DataType outputDataType = qobject_cast( mLayers.at( 0 ) )->dataProvider()->dataType( 1 ); + double outputNoData = 0.0; + if ( qobject_cast( mLayers.at( 0 ) )->dataProvider()->sourceHasNoDataValue( 1 ) ) + { + outputNoData = qobject_cast( mLayers.at( 0 ) )->dataProvider()->sourceNoDataValue( 1 ); + } + else + { + outputNoData = -FLT_MAX; + } + QgsRectangle outputExtent; if ( parameters.value( QStringLiteral( "EXTENT" ) ).isValid() ) { @@ -171,7 +183,7 @@ QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶met 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, 0 ); + provider->setNoDataValue( 1, outputNoData ); const double step = rows > 0 ? 100.0 / rows : 1; for ( int row = 0; row < rows; row++ ) @@ -227,7 +239,7 @@ QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶met } else { - block.setValue( 0, col, 0 ); + block.setValue( 0, col, outputNoData ); } } From e0b02d02e704a437b93717bb8793d83bdcd72008 Mon Sep 17 00:00:00 2001 From: Mathieu Pellerin Date: Fri, 14 Feb 2025 08:54:27 +0700 Subject: [PATCH 3/9] Add a parameter to define NoData handling behavior --- src/analysis/processing/qgsalgorithmrasterrank.cpp | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/analysis/processing/qgsalgorithmrasterrank.cpp b/src/analysis/processing/qgsalgorithmrasterrank.cpp index 191df6d38b99..8c7e4db9c53d 100644 --- a/src/analysis/processing/qgsalgorithmrasterrank.cpp +++ b/src/analysis/processing/qgsalgorithmrasterrank.cpp @@ -60,6 +60,7 @@ void QgsRasterRankAlgorithm::initAlgorithm( const QVariantMap & ) { addParameter( new QgsProcessingParameterMultipleLayers( QStringLiteral( "LAYERS" ), QObject::tr( "Input layers" ), Qgis::ProcessingSourceType::Raster ) ); addParameter( new QgsProcessingParameterNumber( QStringLiteral( "RANK" ), QObject::tr( "Rank" ), Qgis::ProcessingNumberParameterType::Integer, 1 ) ); + 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( 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" ) ); @@ -79,6 +80,12 @@ void QgsRasterRankAlgorithm::initAlgorithm( const QVariantMap & ) bool QgsRasterRankAlgorithm::prepareAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) { + if ( parameterAsInt( parameters, QStringLiteral( "RANK" ), context ) == 0 ) + { + feedback->reportError( QObject::tr( "Rank value must be greater than zero" ), false ); + return false; + } + const QList layers = parameterAsLayerList( parameters, QStringLiteral( "LAYERS" ), context ); for ( const QgsMapLayer *layer : std::as_const( layers ) ) @@ -109,6 +116,7 @@ QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶met } const int rank = parameterAsInt( parameters, QStringLiteral( "RANK" ), context ); + QgsCoordinateReferenceSystem outputCrs; if ( parameters.value( QStringLiteral( "CRS" ) ).isValid() ) { @@ -130,6 +138,7 @@ QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶met { outputNoData = -FLT_MAX; } + const bool outputNoDataOverride = parameterAsInt( parameters, QStringLiteral( "NODATA_HANDLING" ), context ) == 1; QgsRectangle outputExtent; if ( parameters.value( QStringLiteral( "EXTENT" ) ).isValid() ) @@ -230,6 +239,11 @@ QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶met { values << value; } + else if ( outputNoDataOverride ) + { + values.clear(); + break; + } } std::sort( values.begin(), values.end() ); From 72ac0df1f057e4a45f7089268ea946de5cf6f898 Mon Sep 17 00:00:00 2001 From: Mathieu Pellerin Date: Fri, 14 Feb 2025 09:10:18 +0700 Subject: [PATCH 4/9] Allow for negative rank value for reverse ranking --- src/analysis/processing/qgsalgorithmrasterrank.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/analysis/processing/qgsalgorithmrasterrank.cpp b/src/analysis/processing/qgsalgorithmrasterrank.cpp index 8c7e4db9c53d..50737262f2c4 100644 --- a/src/analysis/processing/qgsalgorithmrasterrank.cpp +++ b/src/analysis/processing/qgsalgorithmrasterrank.cpp @@ -247,9 +247,9 @@ QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶met } std::sort( values.begin(), values.end() ); - if ( values.size() >= rank ) + if ( values.size() >= std::abs( rank ) ) { - block.setValue( 0, col, values.at( rank - 1 ) ); + block.setValue( 0, col, values.at( rank > 0 ? rank - 1 : values.size() + rank ) ); } else { From 739c17057f40454ce799d0d4785b289e35fc2bf0 Mon Sep 17 00:00:00 2001 From: Mathieu Pellerin Date: Fri, 14 Feb 2025 09:50:09 +0700 Subject: [PATCH 5/9] Handle non-square cells --- .../processing/qgsalgorithmrasterrank.cpp | 30 ++++++++++++------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/src/analysis/processing/qgsalgorithmrasterrank.cpp b/src/analysis/processing/qgsalgorithmrasterrank.cpp index 50737262f2c4..6657e7a062d4 100644 --- a/src/analysis/processing/qgsalgorithmrasterrank.cpp +++ b/src/analysis/processing/qgsalgorithmrasterrank.cpp @@ -150,7 +150,8 @@ QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶met outputExtent = QgsProcessingUtils::combineLayerExtents( mLayers, outputCrs, context ); } - double minCellSize = 1e9; + double minCellSizeX = 1e9; + double minCellSizeY = 1e9; std::map> inputBlocks; for ( QgsMapLayer *layer : mLayers ) { @@ -163,23 +164,30 @@ QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶met extent = ct.transformBoundingBox( extent ); } - double cellSize = ( extent.xMaximum() - extent.xMinimum() ) / rLayer->width(); - if ( cellSize < minCellSize ) + double cellSizeX = ( extent.xMaximum() - extent.xMinimum() ) / rLayer->width(); + if ( cellSizeX < minCellSizeX ) { - minCellSize = cellSize; + minCellSizeX = cellSizeX; + } + double cellSizeY = ( extent.yMaximum() - extent.yMinimum() ) / rLayer->height(); + if ( cellSizeY < minCellSizeY ) + { + minCellSizeY = cellSizeY; } inputBlocks[rLayer->id()] = std::make_unique(); } - double outputCellSize = parameterAsDouble( parameters, QStringLiteral( "CELL_SIZE" ), context ); - if ( outputCellSize == 0 ) + double outputCellSizeX = parameterAsDouble( parameters, QStringLiteral( "CELL_SIZE" ), context ); + double outputCellSizeY = outputCellSizeX; + if ( outputCellSizeX == 0 ) { - outputCellSize = minCellSize; + outputCellSizeX = minCellSizeX; + outputCellSizeY = minCellSizeY; } - double cols = std::round( ( outputExtent.xMaximum() - outputExtent.xMinimum() ) / outputCellSize ); - double rows = std::round( ( outputExtent.yMaximum() - outputExtent.yMinimum() ) / outputCellSize ); + double cols = std::round( ( outputExtent.xMaximum() - outputExtent.xMinimum() ) / outputCellSizeX ); + double rows = std::round( ( outputExtent.yMaximum() - outputExtent.yMinimum() ) / outputCellSizeY ); const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context ); const QFileInfo fi( outputFile ); @@ -207,8 +215,8 @@ QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶met // Calculates the rect for a single row read QgsRectangle rowExtent( outputExtent ); - rowExtent.setYMaximum( rowExtent.yMaximum() - outputCellSize * row ); - rowExtent.setYMinimum( rowExtent.yMaximum() - outputCellSize ); + rowExtent.setYMaximum( rowExtent.yMaximum() - outputCellSizeY * row ); + rowExtent.setYMinimum( rowExtent.yMaximum() - outputCellSizeY ); for ( QgsMapLayer *layer : mLayers ) { From 816bbe34d713ecff859361330e25201c39dd1faf Mon Sep 17 00:00:00 2001 From: Mathieu Pellerin Date: Fri, 14 Feb 2025 10:26:15 +0700 Subject: [PATCH 6/9] Add test coverage --- .../src/analysis/testqgsprocessingalgspt1.cpp | 123 ++++++++++++++++++ .../rasterRank_testcase1.tif | Bin 0 -> 10504 bytes .../rasterRank_testcase2.tif | Bin 0 -> 10504 bytes .../rasterRank_testcase3.tif | Bin 0 -> 10504 bytes tests/testdata/raster/rank1.tif | Bin 0 -> 10374 bytes tests/testdata/raster/rank2.tif | Bin 0 -> 10374 bytes tests/testdata/raster/rank3.tif | Bin 0 -> 10374 bytes tests/testdata/raster/rank4.tif | Bin 0 -> 10374 bytes 8 files changed, 123 insertions(+) create mode 100644 tests/testdata/control_images/expected_rasterRank/rasterRank_testcase1.tif create mode 100644 tests/testdata/control_images/expected_rasterRank/rasterRank_testcase2.tif create mode 100644 tests/testdata/control_images/expected_rasterRank/rasterRank_testcase3.tif create mode 100644 tests/testdata/raster/rank1.tif create mode 100644 tests/testdata/raster/rank2.tif create mode 100644 tests/testdata/raster/rank3.tif create mode 100644 tests/testdata/raster/rank4.tif diff --git a/tests/src/analysis/testqgsprocessingalgspt1.cpp b/tests/src/analysis/testqgsprocessingalgspt1.cpp index bd14f0317ef2..55f27e2f2ed9 100644 --- a/tests/src/analysis/testqgsprocessingalgspt1.cpp +++ b/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,126 @@ void TestQgsProcessingAlgsPt1::createConstantRaster() } } +void TestQgsProcessingAlgsPt1::rasterRank_data() +{ + QTest::addColumn( "expectedRaster" ); + QTest::addColumn( "rank" ); + QTest::addColumn( "nodataHandling" ); + + /* + * Testcase 1 + * + * inputExtent = from "/raster/band1_int16_noct_epsg4326.tif" + * crs = EPSG:4326 + * pixelSize = 1.0 + * constantValue = 12 + * Byte Raster Layer + * + */ + QTest::newRow( "testcase 1" ) + << QStringLiteral( "/rasterRank_testcase1.tif" ) + << 2 + << 0; + + /* + * Testcase 2 + * + * inputExtent = from "/raster/band1_int16_noct_epsg4326.tif" + * crs = EPSG:4326 + * pixelSize = 1.0 + * constantValue = -1 + * Byte Raster Layer + * + */ + QTest::newRow( "testcase 2" ) + << QStringLiteral( "/rasterRank_testcase2.tif" ) + << -2 + << 0; + + /* + * Testcase 3 + * + * inputExtent = from "/raster/band1_int16_noct_epsg4326.tif" + * crs = EPSG:4326 + * pixelSize = 1.0 + * constantValue = -1 + * Byte Raster Layer + * + */ + QTest::newRow( "testcase 3" ) + << QStringLiteral( "/rasterRank_testcase3.tif" ) + << 2 + << 1; +} + +void TestQgsProcessingAlgsPt1::rasterRank() +{ + QFETCH( QString, expectedRaster ); + QFETCH( int, rank ); + QFETCH( int, nodataHandling ); + + //prepare input params + std::unique_ptr alg( QgsApplication::processingRegistry()->createAlgorithmById( QStringLiteral( "native:rasterrank" ) ) ); + + const QString testDataPath( TEST_DATA_DIR ); //defined in CMakeLists.txt + + QVariantMap parameters; + + parameters.insert( QStringLiteral( "LAYERS" ), QStringList() << testDataPath + "/raster/rank1.tif" << testDataPath + "/raster/rank2.tif" << testDataPath + "/raster/rank3.tif" << testDataPath + "/raster/rank4.tif" ); + parameters.insert( QStringLiteral( "RANK" ), rank ); + 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( testDataPath + "/control_images/expected_rasterRank" + expectedRaster, "expectedDataset", "gdal" ); + std::unique_ptr expectedInterface( expectedRasterLayer->dataProvider()->clone() ); + QgsRasterIterator expectedIter( expectedInterface.get() ); + expectedIter.startRasterRead( 1, expectedRasterLayer->width(), expectedRasterLayer->height(), expectedInterface->extent() ); + + //run algorithm... + auto context = std::make_unique(); + results = alg->run( parameters, *context, &feedback, &ok ); + QVERIFY( ok ); + + //...and check results with expected datasets + auto outputRaster = std::make_unique( results.value( QStringLiteral( "OUTPUT" ) ).toString(), "output", "gdal" ); + std::unique_ptr 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 outputRasterBlock; + std::unique_ptr 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() { diff --git a/tests/testdata/control_images/expected_rasterRank/rasterRank_testcase1.tif b/tests/testdata/control_images/expected_rasterRank/rasterRank_testcase1.tif new file mode 100644 index 0000000000000000000000000000000000000000..3aab8d01faa4e93250796936ccbe7fc71c2a5e4d GIT binary patch literal 10504 zcmeI2ze>bF5XSeCLq${&Q4s_!#Y(Ui6e~*+d;qbqun_D8|MXfZHda@Ol{S`w*w|U= z1Na1%R!%mXWM;CW4s$t+351(`+5LVyUoKZ1yS6s(2Irg`cCMAWmUC_!&K9bZP!E8f zhB_P9hw9N*+a;h~!#R%hiBPYDJ_&U;COf9lmmb$PqGGPVEuaelYB5*uThC143y0+xU!Ujc{Vhz!qfwA+zef0YU==2gl}|EqmAzVaew-^daIVN|goLM6)lSs&n>`w#?aPsX8s zxgphLzIYb(^NfguX9((=BNG&|!e1=!UKa!lXH)pH`869Zy5ycI6GWzL60=MjVK4|= zlI)T+lEOg*jjt8^HK{RV>*lk_vNHy0gpU?pV@LxGkHCEhmtCJ$HPQehO~O%7>+2CK zxg|7FFp|VeV#CP9%+SYFF+3zPR4xSs>H}yxuV1yO{9#i@byuYd2+An_PE;Ylznif> z5>-fnV0~ypN<>c>YH3w%;FmCeRg|80`M@(4dI>clxQ*U0be;(5qf3o6z@Uik^uZWi zf%Jmm_$#9D8nQlIBMc5v@539?w6qz)k3f_er5%QpqP#%Q9GM^pQl(|dU5|L9 zh)vBgu{)+&0+xU!U9f0?+m++?@fSsA${-1ULK7dcL1Mm&% V!%pIN4=#kTKH>9f);X=+fM0$GLXQ9d literal 0 HcmV?d00001 diff --git a/tests/testdata/control_images/expected_rasterRank/rasterRank_testcase2.tif b/tests/testdata/control_images/expected_rasterRank/rasterRank_testcase2.tif new file mode 100644 index 0000000000000000000000000000000000000000..d763018c6097e1c599d4acd9de73e09b462de420 GIT binary patch literal 10504 zcmeI2Jxjzu5QaA!IYrI`5fwquQmh1P;SMWH5&QvSVPPTI3x4!kDK=JDiPhO$DTvK= zR{8_{1(sIMCYxktHyBMOr@|~GnaqAX_r9Bix3Ep4{EUq^AwLVTX=BPfYO{sn)@f6PP~U^?emc5n9lpF@-miA@+DJ|u5C_BoaX=gp2gCt!KpgnC z17TEhng_PzRQ)s+Ghk0UCWOGekXI}zf!7>*DQ9YsQhz+gf6I*3_33?G88tRXO% z5dN)>bkxSsI-#&3JRurNL}5g@vu+p>gp-jA(JYTc^fT!bJ`$4r(}!AB3biH>YgVNp zz0NQwh8TvKbgx~7MD)8eiY2f@Vg&51mh|>u6p}J*5ZxQ`u7WwFtBc5F@Xg3q3?{cc z1lCa(15F3lrRC}PpNs+-EE7`1;DERcIvYHb?lS6lZ+95Ae1CWhR|t9Tlsln}Xbd_f z<%T7_2*cJ82x(pJ1T)eRFeFWu4-f+lf*>ecHcT1L04c@EmpQeJysj}OuSRm>fH)uy zhy&t)I3Ny)1LDBnaG=Gq0fi=hO$BG*pS$Qj4TXHK(SRoEQ%y9 z24x2V#l?YakQxw>K~fWp#O7<}VPFQzZvx_`79IvRAbkpmk-?5r8Bc`6raSbnIehVG z-2?|@Ic)4D1-3~uH#s{>7GM7~ZL_f>HZ`MSqaiRF0;3@?8UmvsFd71*Aut*OqaiRF z0;3@?8UmvsFd71*Aut*OgE0h{z<{%fXa(Q_pgFJ=fH#2b#wASP)q@7mp!x-*=P58SP-`_H0Pb8fhX4Qo literal 0 HcmV?d00001 diff --git a/tests/testdata/raster/rank1.tif b/tests/testdata/raster/rank1.tif new file mode 100644 index 0000000000000000000000000000000000000000..648f83d811a4e55a54de1507ee831242675a1d2c GIT binary patch literal 10374 zcmeI1ze~eF6vy8sq0LaNBK}0Q(n%ex6o(E%OBV|+B6M^R2L}g1a4Ld7Ae%Z>+&c{Q0;^rTMWEK) z+D_?iZRaw9dPEqY;oK|UV}bc`x8B7039Ri*XKp_`6Ub}@>Y=k7=$LuO`}8~q863ax zzO!^F;5dumuBE{w*4tR0bTp@D37_rj2zgXZa>02G{@_U;}J`4X^=4284Gh>NnWN132U0YZWs+2z*e{ zv$sb-;SV4o2<0aUWvJB^j!?>oKq$2p$xwHaXhbRuRgjotDEmZW8EOs^k6=dd`%vOa zg)A9@8U6>tu6!4=ie1SJ`S014??PJ1WJwzk0{?(P*rhxoIfuLfB@qvybdYCa8S>?% zL>dhIM^5BL+9V|uNl^A@BJLs*LfIiPlLjN|AjB;aW1%M@t|iApAmp{=-~fcYGC{Bp ml3-<$Xn>MrS<>t;Ch029vjH~12G{@_U;}J`4X}aq4Ez8w>nwx- literal 0 HcmV?d00001 diff --git a/tests/testdata/raster/rank2.tif b/tests/testdata/raster/rank2.tif new file mode 100644 index 0000000000000000000000000000000000000000..ca771316b83d756d90d493da1da7332c279adb74 GIT binary patch literal 10374 zcmeHIJxc>Y5S_b(;7XAX~j=Q+Km)qM7d!jNIHZ!yH=Do~IAocnbH~`QMU}eCg~?io)XbKf82D|kPExrb?tZRbZA%qo+29Ty;jm{&YUCk4o1{J?YL zW}t*|6os3v4TG4kV}97MeaxY-egY%l-l)N=|<~Ocx z7Z-+h-##zyDhnah@*T+oWYH7v6+}DRbc+bGoTlqQ5prmb1*a2G@S@{n{Wv+XxHMZW zSEfrdrRBBi+pf`}-Qf=U|+)0?!DFl|?R4WpaHNSNFtR>I&S@fyZNBoaTOVMs)_hN*8| zu3=`8f1;7Fh7qrlG)#z)hWbAxc|~4Dego>3#3Cb{;zjJFdPAg))aq2R$Oxx+5j&~g z5Glo4ok@ZeNrhT6sooMe#ZsL}f)q*lS~959;z@-{he(1ISou;ir_?6?>Iq=6F_yV8IF=GG# literal 0 HcmV?d00001 diff --git a/tests/testdata/raster/rank3.tif b/tests/testdata/raster/rank3.tif new file mode 100644 index 0000000000000000000000000000000000000000..4e1409f966988bc3a67de37f15b2befd6e9f83ce GIT binary patch literal 10374 zcmeH|Jxc>Y5QgVINUjJ*5kC=)q*4nb!6Jo_kj8+G2q`VZ!oorjtb*VNT&GrweM&*1 zjlI}d2sZi$tS$Tn&Yc%HlAGPzV@(j}!e%Br@62=iZjoAT9vlD|1TZWxFx$pugmN>; zJ(U&Yo-P&2hYTA`+spY*IE~?yYQH z-&HFU`|n?u_r*$EYWR+%0c`Zd`vlRnTW%IbhSPE#$U+LuG2ko@8814H?Z>%%b#JHQSoIbi)fN2&LglEGUyid3I$Q0BJHQUG1MC1h&_4&h0ZwiY5QgVINUjJ*5kC=)q*4nb!6Jo_kj8+G2q`VZ!oorjtb*VNT&GqM`;>x2 z8+);_5Nz}hSX=lDoIUQ~NV2=PVegP*E^KCQ_uYMV-d(sxV;LL(7zHpaFfiN4W`uIn z&zZ{dbEeCL@-f2(6Z_t=O@#U<{d^nyr!Z%l&cb1KE~MEB8Yxj_Ivg!guzZVBhQ_!1LR|_jvQRdVT-$@p^hw^#%ohouiGt zyT^KM=J@Nk{Zy_6Qp0m33t*!s?B_?*?z(vt8BW)AAP*Td$AGIMq& z@p^F)i7>Q5_+{uj5fMXEAwooI7+MMuTczjJi4pxfCi6i*CxRbo5kOo(uk-KrG-Rd0Yo%S zX^hWMT1G79LBzAvMw4K0u9n(HrRqVINGGWZ`XQ<-M9}v_B1{=V5~d6>J>ryr5>ZYv ja)2Bl2gm_(fE*wP$N_SI93ThC0djyGAP2~SbUW||qN+PF literal 0 HcmV?d00001 From f3b0e31f7bbe5f53cb26afc5c53b847620a25ff3 Mon Sep 17 00:00:00 2001 From: Mathieu Pellerin Date: Fri, 14 Feb 2025 16:28:44 +0700 Subject: [PATCH 7/9] Allow for multiple ranks to be computed at once, with the output raster adopting one band per rank --- .../processing/qgsalgorithmrasterrank.cpp | 70 ++++++++++++------- .../processing/qgsalgorithmrasterrank.h | 1 + .../src/analysis/testqgsprocessingalgspt1.cpp | 41 ++++------- 3 files changed, 58 insertions(+), 54 deletions(-) diff --git a/src/analysis/processing/qgsalgorithmrasterrank.cpp b/src/analysis/processing/qgsalgorithmrasterrank.cpp index 6657e7a062d4..0b39a8244e3f 100644 --- a/src/analysis/processing/qgsalgorithmrasterrank.cpp +++ b/src/analysis/processing/qgsalgorithmrasterrank.cpp @@ -59,7 +59,7 @@ QgsRasterRankAlgorithm *QgsRasterRankAlgorithm::createInstance() const void QgsRasterRankAlgorithm::initAlgorithm( const QVariantMap & ) { addParameter( new QgsProcessingParameterMultipleLayers( QStringLiteral( "LAYERS" ), QObject::tr( "Input layers" ), Qgis::ProcessingSourceType::Raster ) ); - addParameter( new QgsProcessingParameterNumber( QStringLiteral( "RANK" ), QObject::tr( "Rank" ), Qgis::ProcessingNumberParameterType::Integer, 1 ) ); + addParameter( new QgsProcessingParameterString( QStringLiteral( "RANKS" ), QObject::tr( "Rank (separate multiple ranks using commas)" ), 1 ) ); 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( QStringLiteral( "EXTENT" ), QObject::tr( "Output extent" ), QVariant(), true ); @@ -80,14 +80,24 @@ void QgsRasterRankAlgorithm::initAlgorithm( const QVariantMap & ) bool QgsRasterRankAlgorithm::prepareAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) { - if ( parameterAsInt( parameters, QStringLiteral( "RANK" ), context ) == 0 ) + const QStringList rankStrings = parameterAsString( parameters, QStringLiteral( "RANKS" ), context ).split( QLatin1String( "," ) ); + for ( const QString &rankString : rankStrings ) { - feedback->reportError( QObject::tr( "Rank value must be greater than zero" ), false ); + bool ok = false; + const int rank = rankString.toInt( &ok ); + if ( ok && rank != 0 ) + { + mRanks << rank; + } + } + + if ( mRanks.isEmpty() ) + { + feedback->reportError( QObject::tr( "No valid non-zero rank value(s) provided" ), false ); return false; } const QList layers = parameterAsLayerList( parameters, QStringLiteral( "LAYERS" ), context ); - for ( const QgsMapLayer *layer : std::as_const( layers ) ) { if ( !qobject_cast( layer ) || !layer->dataProvider() ) @@ -115,8 +125,6 @@ QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶met layer->moveToThread( QThread::currentThread() ); } - const int rank = parameterAsInt( parameters, QStringLiteral( "RANK" ), context ); - QgsCoordinateReferenceSystem outputCrs; if ( parameters.value( QStringLiteral( "CRS" ) ).isValid() ) { @@ -164,16 +172,8 @@ QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶met extent = ct.transformBoundingBox( extent ); } - double cellSizeX = ( extent.xMaximum() - extent.xMinimum() ) / rLayer->width(); - if ( cellSizeX < minCellSizeX ) - { - minCellSizeX = cellSizeX; - } - double cellSizeY = ( extent.yMaximum() - extent.yMinimum() ) / rLayer->height(); - if ( cellSizeY < minCellSizeY ) - { - minCellSizeY = cellSizeY; - } + minCellSizeX = std::min( minCellSizeX, ( extent.xMaximum() - extent.xMinimum() ) / rLayer->width() ); + minCellSizeY = std::min( minCellSizeY, ( extent.yMaximum() - extent.yMinimum() ) / rLayer->height() ); inputBlocks[rLayer->id()] = std::make_unique(); } @@ -195,13 +195,20 @@ QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶met auto writer = std::make_unique( outputFile ); writer->setOutputFormat( outputFormat ); - std::unique_ptr provider( writer->createOneBandRaster( outputDataType, cols, rows, outputExtent, outputCrs ) ); + std::unique_ptr 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::vector> outputBlocks; + for ( int i = 0; i < mRanks.size(); i++ ) + { + outputBlocks.push_back( std::make_unique() ); + } + const double step = rows > 0 ? 100.0 / rows : 1; for ( int row = 0; row < rows; row++ ) { @@ -210,8 +217,11 @@ QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶met break; } - QgsRasterBlock block( outputDataType, cols, 1 ); - block.setNoDataValue( 0 ); + for ( int i = 0; i < mRanks.size(); i++ ) + { + outputBlocks[i].reset( new QgsRasterBlock( outputDataType, cols, 1 ) ); + outputBlocks[i]->setNoDataValue( outputNoData ); + } // Calculates the rect for a single row read QgsRectangle rowExtent( outputExtent ); @@ -253,19 +263,25 @@ QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶met break; } } - std::sort( values.begin(), values.end() ); - if ( values.size() >= std::abs( rank ) ) - { - block.setValue( 0, col, values.at( rank > 0 ? rank - 1 : values.size() + rank ) ); - } - else + + for ( int i = 0; i < mRanks.size(); i++ ) { - block.setValue( 0, col, outputNoData ); + if ( values.size() >= std::abs( mRanks[i] ) ) + { + outputBlocks[i]->setValue( 0, col, values.at( mRanks[i] > 0 ? mRanks[i] - 1 : values.size() + mRanks[i] ) ); + } + else + { + outputBlocks[i]->setValue( 0, col, outputNoData ); + } } } - provider->writeBlock( &block, 1, 0, row ); + for ( int i = 0; i < mRanks.size(); i++ ) + { + provider->writeBlock( outputBlocks[i].get(), i + 1, 0, row ); + } feedback->setProgress( row * step ); } diff --git a/src/analysis/processing/qgsalgorithmrasterrank.h b/src/analysis/processing/qgsalgorithmrasterrank.h index 6b298ace9516..895e61d95178 100644 --- a/src/analysis/processing/qgsalgorithmrasterrank.h +++ b/src/analysis/processing/qgsalgorithmrasterrank.h @@ -46,6 +46,7 @@ class QgsRasterRankAlgorithm : public QgsProcessingAlgorithm bool prepareAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override; QVariantMap processAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override; + QList mRanks; QList mLayers; }; diff --git a/tests/src/analysis/testqgsprocessingalgspt1.cpp b/tests/src/analysis/testqgsprocessingalgspt1.cpp index 55f27e2f2ed9..49ec1069c281 100644 --- a/tests/src/analysis/testqgsprocessingalgspt1.cpp +++ b/tests/src/analysis/testqgsprocessingalgspt1.cpp @@ -1903,59 +1903,46 @@ void TestQgsProcessingAlgsPt1::createConstantRaster() void TestQgsProcessingAlgsPt1::rasterRank_data() { QTest::addColumn( "expectedRaster" ); - QTest::addColumn( "rank" ); + QTest::addColumn( "ranks" ); QTest::addColumn( "nodataHandling" ); /* * Testcase 1 - * - * inputExtent = from "/raster/band1_int16_noct_epsg4326.tif" - * crs = EPSG:4326 - * pixelSize = 1.0 - * constantValue = 12 - * Byte Raster Layer - * */ QTest::newRow( "testcase 1" ) << QStringLiteral( "/rasterRank_testcase1.tif" ) - << 2 + << QString::number( 2 ) << 0; /* * Testcase 2 - * - * inputExtent = from "/raster/band1_int16_noct_epsg4326.tif" - * crs = EPSG:4326 - * pixelSize = 1.0 - * constantValue = -1 - * Byte Raster Layer - * */ QTest::newRow( "testcase 2" ) << QStringLiteral( "/rasterRank_testcase2.tif" ) - << -2 + << QString::number( -2 ) << 0; /* * Testcase 3 - * - * inputExtent = from "/raster/band1_int16_noct_epsg4326.tif" - * crs = EPSG:4326 - * pixelSize = 1.0 - * constantValue = -1 - * Byte Raster Layer - * */ QTest::newRow( "testcase 3" ) << QStringLiteral( "/rasterRank_testcase3.tif" ) - << 2 + << 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( int, rank ); + QFETCH( QString, ranks ); QFETCH( int, nodataHandling ); //prepare input params @@ -1966,7 +1953,7 @@ void TestQgsProcessingAlgsPt1::rasterRank() QVariantMap parameters; parameters.insert( QStringLiteral( "LAYERS" ), QStringList() << testDataPath + "/raster/rank1.tif" << testDataPath + "/raster/rank2.tif" << testDataPath + "/raster/rank3.tif" << testDataPath + "/raster/rank4.tif" ); - parameters.insert( QStringLiteral( "RANK" ), rank ); + parameters.insert( QStringLiteral( "RANKS" ), ranks ); parameters.insert( QStringLiteral( "NODATA_HANDLING" ), nodataHandling ); parameters.insert( QStringLiteral( "OUTPUT" ), QgsProcessing::TEMPORARY_OUTPUT ); From 85f5c813249eff0cfb565b1f37258089d512affc Mon Sep 17 00:00:00 2001 From: Mathieu Pellerin Date: Sat, 15 Feb 2025 09:11:55 +0700 Subject: [PATCH 8/9] Add missing test reference image --- .../rasterRank_testcase4.tif | Bin 0 -> 20772 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 tests/testdata/control_images/expected_rasterRank/rasterRank_testcase4.tif diff --git a/tests/testdata/control_images/expected_rasterRank/rasterRank_testcase4.tif b/tests/testdata/control_images/expected_rasterRank/rasterRank_testcase4.tif new file mode 100644 index 0000000000000000000000000000000000000000..ddf528e368f5eb406ba4eda4c8c67c912511016f GIT binary patch literal 20772 zcmeI4ziSjh6vyArEaw)7XApmZkb^W)3nRfIg>WH_0l{C8KR_%}iZ&KO{1LcLQz?Q) zij;z!HkN|ehzNpUYvEsDW$ByUoy^Sc%6Vb<5r!tPao zBLE491ORU1m1s6Yo1rhmlb3$cuBX{^NeUc&-|;$D`bD4PrM-yu3A9CGxb{B?BW#pv1V zkKaB$*cn+HY0kYnmu|g$yS94j-jCnU-z~4oTOuNWLjY+J9RlZ3vUD;Y;GHC!j58QO zFBAT!3(yy3mZf=lW?}8}g_Wh{bHnq)D;p~px61049G*&2fD|AFNC8rS6d(mi0aAbz zAO%PPQh*d71xNu>fD|AFT2LV6>YQ%6G%?TVdL48tMIkNUj3z-*^R;j5o{ZfrU9;OL zQs;N$u2BxcZ^p0JsG zmAB1r23^$T+bJ)NOh_8rNFjYT3gOd^PIJVy+FjXb&d1DgSf8rc-dgT_>}P~HT6=x@ zvQg&4lZ|**h@T@T@p8mM(#X+3F)JJZ^3&T@wzdSk6Z!KYy1L?*e(&(D-E8J)B55sD z??gTh)lq71D^;tImYc(j%0^K5%yyw_E{=LUE|Pl75FR^+JrX;Ky=~jxI)^`kI*Y%3 zTX#*4X=l}GrtL%RNC8rS6d(mi0aAbzAO%PPQh*d71xNu>fD|AFNC8rS6d(nZ0&}>H z-2d;(cmG=RzI=DF&A;r+w~g)PzI?CHmizL3L0j(2*S`+%XpiN4`Yn}EqxsSN?%*N* PpW25Zx<&x8@6q=cnh Date: Tue, 25 Feb 2025 11:19:21 +0700 Subject: [PATCH 9/9] Optimize further by iterating using QgsRasterIterator --- .../processing/qgsalgorithmrasterrank.cpp | 117 ++++++++++-------- 1 file changed, 65 insertions(+), 52 deletions(-) diff --git a/src/analysis/processing/qgsalgorithmrasterrank.cpp b/src/analysis/processing/qgsalgorithmrasterrank.cpp index 0b39a8244e3f..ab22320ae448 100644 --- a/src/analysis/processing/qgsalgorithmrasterrank.cpp +++ b/src/analysis/processing/qgsalgorithmrasterrank.cpp @@ -160,7 +160,6 @@ QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶met double minCellSizeX = 1e9; double minCellSizeY = 1e9; - std::map> inputBlocks; for ( QgsMapLayer *layer : mLayers ) { QgsRasterLayer *rLayer = qobject_cast( layer ); @@ -174,8 +173,6 @@ QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶met minCellSizeX = std::min( minCellSizeX, ( extent.xMaximum() - extent.xMinimum() ) / rLayer->width() ); minCellSizeY = std::min( minCellSizeY, ( extent.yMaximum() - extent.yMinimum() ) / rLayer->height() ); - - inputBlocks[rLayer->id()] = std::make_unique(); } double outputCellSizeX = parameterAsDouble( parameters, QStringLiteral( "CELL_SIZE" ), context ); @@ -202,6 +199,29 @@ QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶met throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( outputFile, provider->error().message( QgsErrorMessage::Text ) ) ); provider->setNoDataValue( 1, outputNoData ); + std::map> inputInterfaces; + std::map> inputBlocks; + for ( QgsMapLayer *layer : mLayers ) + { + QgsRasterLayer *rLayer = qobject_cast( layer ); + if ( rLayer->crs() != outputCrs ) + { + QgsRasterProjector *projector = new QgsRasterProjector(); + projector->setCrs( rLayer->crs(), outputCrs, context.transformContext() ); + projector->setInput( rLayer->dataProvider() ); + projector->setPrecision( QgsRasterProjector::Exact ); + inputInterfaces[rLayer->id()].reset( projector ); + } + else + { + inputInterfaces[rLayer->id()].reset( rLayer->dataProvider()->clone() ); + } + inputBlocks[rLayer->id()] = std::make_unique(); + } + + std::unique_ptr rasterIterator = std::make_unique( inputInterfaces[mLayers.at( 0 )->id()].get() ); + rasterIterator->startRasterRead( 1, cols, rows, outputExtent ); + int rasterIteratorCount = static_cast( rasterIterator->blockCount() ); std::vector> outputBlocks; for ( int i = 0; i < mRanks.size(); i++ ) @@ -209,80 +229,73 @@ QVariantMap QgsRasterRankAlgorithm::processAlgorithm( const QVariantMap ¶met outputBlocks.push_back( std::make_unique() ); } - const double step = rows > 0 ? 100.0 / rows : 1; - for ( int row = 0; row < rows; row++ ) + + const double step = rasterIteratorCount > 0 ? 100.0 / rasterIteratorCount : 0; + for ( int i = 0; i < rasterIteratorCount; i++ ) { if ( feedback->isCanceled() ) { break; } + feedback->setProgress( i * step ); - for ( int i = 0; i < mRanks.size(); i++ ) + 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 ) { - outputBlocks[i].reset( new QgsRasterBlock( outputDataType, cols, 1 ) ); - outputBlocks[i]->setNoDataValue( outputNoData ); + inputBlocks[inputInterface.first].reset( inputInterface.second->block( 1, blockExtent, iterCols, iterRows ) ); } - // Calculates the rect for a single row read - QgsRectangle rowExtent( outputExtent ); - rowExtent.setYMaximum( rowExtent.yMaximum() - outputCellSizeY * row ); - rowExtent.setYMinimum( rowExtent.yMaximum() - outputCellSizeY ); - - for ( QgsMapLayer *layer : mLayers ) + for ( int i = 0; i < mRanks.size(); i++ ) { - QgsRasterLayer *rLayer = qobject_cast( layer ); - - if ( rLayer->crs() != outputCrs ) - { - QgsRasterProjector proj; - proj.setCrs( rLayer->crs(), outputCrs, context.transformContext() ); - proj.setInput( rLayer->dataProvider() ); - proj.setPrecision( QgsRasterProjector::Exact ); - inputBlocks[rLayer->id()].reset( proj.block( 1, rowExtent, cols, 1 ) ); - } - else - { - inputBlocks[rLayer->id()].reset( rLayer->dataProvider()->block( 1, rowExtent, cols, 1 ) ); - } + outputBlocks[i].reset( new QgsRasterBlock( outputDataType, iterCols, iterRows ) ); + outputBlocks[i]->setNoDataValue( outputNoData ); } - for ( int col = 0; col < cols; col++ ) + for ( int row = 0; row < iterRows; row++ ) { - QList values; - for ( const auto &inputBlock : inputBlocks ) + for ( int col = 0; col < iterCols; col++ ) { - bool isNoData = false; - const double value = inputBlock.second->valueAndNoData( 0, col, isNoData ); - if ( !isNoData ) - { - values << value; - } - else if ( outputNoDataOverride ) + QList values; + for ( const auto &inputBlock : inputBlocks ) { - values.clear(); - break; + bool isNoData = false; + const double value = inputBlock.second->valueAndNoData( row, col, isNoData ); + if ( !isNoData ) + { + values << value; + } + else if ( outputNoDataOverride ) + { + values.clear(); + break; + } } - } - std::sort( values.begin(), values.end() ); + std::sort( values.begin(), values.end() ); - for ( int i = 0; i < mRanks.size(); i++ ) - { - if ( values.size() >= std::abs( mRanks[i] ) ) - { - outputBlocks[i]->setValue( 0, col, values.at( mRanks[i] > 0 ? mRanks[i] - 1 : values.size() + mRanks[i] ) ); - } - else + for ( int i = 0; i < mRanks.size(); i++ ) { - outputBlocks[i]->setValue( 0, col, outputNoData ); + if ( values.size() >= std::abs( mRanks[i] ) ) + { + outputBlocks[i]->setValue( row, col, values.at( mRanks[i] > 0 ? mRanks[i] - 1 : values.size() + mRanks[i] ) ); + } + else + { + outputBlocks[i]->setValue( row, col, outputNoData ); + } } } } for ( int i = 0; i < mRanks.size(); i++ ) { - provider->writeBlock( outputBlocks[i].get(), i + 1, 0, row ); + provider->writeBlock( outputBlocks[i].get(), i + 1, iterLeft, iterTop ); } - feedback->setProgress( row * step ); } qDeleteAll( mLayers );