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..ab22320ae448 --- /dev/null +++ b/src/analysis/processing/qgsalgorithmrasterrank.cpp @@ -0,0 +1,309 @@ +/*************************************************************************** + 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 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 ); + 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 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; + } + } + + 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() ) + 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() ); + } + + 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 ); + 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; + } + 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( mLayers, outputCrs, context ); + } + + double minCellSizeX = 1e9; + double minCellSizeY = 1e9; + 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 ); + } + + minCellSizeX = std::min( minCellSizeX, ( extent.xMaximum() - extent.xMinimum() ) / rLayer->width() ); + minCellSizeY = std::min( minCellSizeY, ( extent.yMaximum() - extent.yMinimum() ) / rLayer->height() ); + } + + double outputCellSizeX = parameterAsDouble( parameters, QStringLiteral( "CELL_SIZE" ), context ); + double outputCellSizeY = outputCellSizeX; + if ( outputCellSizeX == 0 ) + { + outputCellSizeX = minCellSizeX; + outputCellSizeY = minCellSizeY; + } + + 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 ); + const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() ); + + auto writer = std::make_unique( outputFile ); + writer->setOutputFormat( outputFormat ); + 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::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++ ) + { + outputBlocks.push_back( std::make_unique() ); + } + + + const double step = rasterIteratorCount > 0 ? 100.0 / rasterIteratorCount : 0; + for ( int i = 0; i < rasterIteratorCount; i++ ) + { + if ( feedback->isCanceled() ) + { + break; + } + feedback->setProgress( i * 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++ ) + { + QList values; + for ( const auto &inputBlock : inputBlocks ) + { + 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() ); + + for ( int i = 0; i < mRanks.size(); i++ ) + { + 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, iterLeft, iterTop ); + } + } + + 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..895e61d95178 --- /dev/null +++ b/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 ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override; + QVariantMap processAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override; + + QList mRanks; + 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() ); diff --git a/tests/src/analysis/testqgsprocessingalgspt1.cpp b/tests/src/analysis/testqgsprocessingalgspt1.cpp index bd14f0317ef2..49ec1069c281 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,113 @@ void TestQgsProcessingAlgsPt1::createConstantRaster() } } +void TestQgsProcessingAlgsPt1::rasterRank_data() +{ + QTest::addColumn( "expectedRaster" ); + QTest::addColumn( "ranks" ); + QTest::addColumn( "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 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( "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( 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 000000000000..3aab8d01faa4 Binary files /dev/null and b/tests/testdata/control_images/expected_rasterRank/rasterRank_testcase1.tif differ 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 000000000000..d763018c6097 Binary files /dev/null and b/tests/testdata/control_images/expected_rasterRank/rasterRank_testcase2.tif differ diff --git a/tests/testdata/control_images/expected_rasterRank/rasterRank_testcase3.tif b/tests/testdata/control_images/expected_rasterRank/rasterRank_testcase3.tif new file mode 100644 index 000000000000..2d14c96b71cd Binary files /dev/null and b/tests/testdata/control_images/expected_rasterRank/rasterRank_testcase3.tif differ 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 000000000000..ddf528e368f5 Binary files /dev/null and b/tests/testdata/control_images/expected_rasterRank/rasterRank_testcase4.tif differ diff --git a/tests/testdata/raster/rank1.tif b/tests/testdata/raster/rank1.tif new file mode 100644 index 000000000000..648f83d811a4 Binary files /dev/null and b/tests/testdata/raster/rank1.tif differ diff --git a/tests/testdata/raster/rank2.tif b/tests/testdata/raster/rank2.tif new file mode 100644 index 000000000000..ca771316b83d Binary files /dev/null and b/tests/testdata/raster/rank2.tif differ diff --git a/tests/testdata/raster/rank3.tif b/tests/testdata/raster/rank3.tif new file mode 100644 index 000000000000..4e1409f96698 Binary files /dev/null and b/tests/testdata/raster/rank3.tif differ diff --git a/tests/testdata/raster/rank4.tif b/tests/testdata/raster/rank4.tif new file mode 100644 index 000000000000..99794480ff80 Binary files /dev/null and b/tests/testdata/raster/rank4.tif differ