From 2716ecd409081f1b1fd76e6152daf12dd7fe7af4 Mon Sep 17 00:00:00 2001 From: Mathieu Pellerin Date: Tue, 25 Feb 2025 11:19:21 +0700 Subject: [PATCH] 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 );