Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
finalize raster stack position algorithms
  • Loading branch information
root676 authored and nyalldawson committed Aug 5, 2020
1 parent 06f9840 commit 58629cb
Show file tree
Hide file tree
Showing 2 changed files with 131 additions and 15 deletions.
138 changes: 127 additions & 11 deletions src/analysis/processing/qgsalgorithmrasterstackposition.cpp
Expand Up @@ -164,18 +164,17 @@ QVariantMap QgsRasterStackPositionAlgorithmBase::processAlgorithm( const QVarian
for ( int col = 0; col < iterCols; col++ )
{
bool noDataInStack = false;
//rewrite getCellValuesFromBlock
std::vector<double> cellValues = QgsRasterAnalysisUtils::getCellValuesFromBlockStack( inputBlocks, row, col, noDataInStack );
int position = findPosition( inputBlocks, row, col, noDataInStack );

if ( noDataInStack && !mIgnoreNoData )
if ( position == -1 || ( noDataInStack && !mIgnoreNoData ) )
{
//output cell will always be NoData if NoData occurs the cellValueStack and NoData is not ignored
//output cell will always be NoData if NoData occurs the current raster cell
//of the input blocks and NoData is not ignored
//this saves unnecessary iterations on the cellValueStack
outputBlock->setValue( row, col, mNoDataValue );
}
else
{
int position = getPosition( cellValues );
outputBlock->setValue( row, col, position );
}
}
Expand Down Expand Up @@ -215,7 +214,7 @@ QStringList QgsRasterStackLowestPositionAlgorithm::tags() const

QString QgsRasterStackLowestPositionAlgorithm::shortHelpString() const
{
return QObject::tr( "The Lowest position algorithm evaluates on a cell-by-cell basis the position "
return QObject::tr( "The lowest position algorithm evaluates on a cell-by-cell basis the position "
"of the raster with the highest value in a stack of rasters. Position counts start "
"with 1 and range to the total number input rasters. The order of the input "
"rasters is relevant for the algorithm. If multiple rasters feature the lowest value, "
Expand All @@ -233,10 +232,68 @@ QgsRasterStackLowestPositionAlgorithm *QgsRasterStackLowestPositionAlgorithm::cr
return new QgsRasterStackLowestPositionAlgorithm();
}

int QgsRasterStackLowestPositionAlgorithm::getPosition( std::vector<QgsRasterBlock> inputBlocks )
int QgsRasterStackLowestPositionAlgorithm::findPosition( std::vector< std::unique_ptr<QgsRasterBlock> > &inputBlocks, int &row, int &col, bool &noDataInRasterBlockStack )
{
int lowestPosition = 0;

//auxiliary variables
int inputBlocksCount = inputBlocks.size();
int currentPosition = 0;
int noDataCount = 0;
double firstValue = mNoDataValue;
bool firstValueIsNoData = true;

while( firstValueIsNoData && ( currentPosition < inputBlocksCount ) )
{
//check if all blocks are nodata/invalid
std::unique_ptr<QgsRasterBlock> &firstBlock = inputBlocks.at(currentPosition);
firstValue = firstBlock->valueAndNoData(row, col, firstValueIsNoData);

if( !firstBlock->isValid() || firstValueIsNoData )
{
noDataInRasterBlockStack = true;
noDataCount++;
}
else
{
lowestPosition = currentPosition;
}
currentPosition++;
}

if( noDataCount == inputBlocksCount )
{
noDataInRasterBlockStack = true;
return -1; //all blocks are NoData
}
else
{
//scan for the lowest value
while( currentPosition < inputBlocksCount)
{
std::unique_ptr< QgsRasterBlock > &currentBlock = inputBlocks.at(currentPosition);

bool currentValueIsNoData = false;
double currentValue = currentBlock->valueAndNoData(row, col, currentValueIsNoData);

if(!currentBlock->isValid() || currentValueIsNoData)
{
noDataInRasterBlockStack = true;
noDataCount++;
}
else
{
if(currentValue < firstValue)
{
firstValue = currentValue;
lowestPosition = currentPosition;
}
}
currentPosition++;
}
}
//the ArcGIS implementation uses 1 for first position value instead of 0 as in standard c++
return static_cast<int>( std::min_element( cellValueStack.begin(), cellValueStack.end() ) - cellValueStack.begin() ) + 1; //therefore... index + 1
return ++lowestPosition; //therefore ++
}

//
Expand All @@ -260,7 +317,7 @@ QStringList QgsRasterStackHighestPositionAlgorithm::tags() const

QString QgsRasterStackHighestPositionAlgorithm::shortHelpString() const
{
return QObject::tr( "The Highest position algorithm evaluates on a cell-by-cell basis the position "
return QObject::tr( "The highest position algorithm evaluates on a cell-by-cell basis the position "
"of the raster with the highest value in a stack of rasters. Position counts start "
"with 1 and range to the total number input rasters. The order of the input "
"rasters is relevant for the algorithm. If multiple rasters feature the highest value, "
Expand All @@ -278,10 +335,69 @@ QgsRasterStackHighestPositionAlgorithm *QgsRasterStackHighestPositionAlgorithm::
return new QgsRasterStackHighestPositionAlgorithm();
}

int QgsRasterStackHighestPositionAlgorithm::getPosition( std::vector<QgsRasterBlock> inputBlocks )
int QgsRasterStackHighestPositionAlgorithm::findPosition( std::vector< std::unique_ptr< QgsRasterBlock> > &inputBlocks, int &row, int &col, bool &noDataInRasterBlockStack )
{
int highestPosition = 0;

//auxiliary variables
int inputBlocksCount = inputBlocks.size();
int currentPosition = 0;
int noDataCount = 0;
double firstValue = mNoDataValue;
bool firstValueIsNoData = true;

while( firstValueIsNoData && ( currentPosition < inputBlocksCount ) )
{
//check if all blocks are nodata/invalid
std::unique_ptr<QgsRasterBlock> &firstBlock = inputBlocks.at(currentPosition);
firstValue = firstBlock->valueAndNoData(row, col, firstValueIsNoData);

if( !firstBlock->isValid() || firstValueIsNoData )
{
noDataInRasterBlockStack = true;
noDataCount++;
}
else
{
highestPosition = currentPosition;
}

currentPosition++;
}

if( noDataCount == inputBlocksCount )
{
noDataInRasterBlockStack = true;
return -1; //all blocks are NoData
}
else
{
//scan for the lowest value
while( currentPosition < inputBlocksCount)
{
std::unique_ptr< QgsRasterBlock > &currentBlock = inputBlocks.at(currentPosition);

bool currentValueIsNoData = false;
double currentValue = currentBlock->valueAndNoData(row, col, currentValueIsNoData);

if(!currentBlock->isValid() || currentValueIsNoData)
{
noDataInRasterBlockStack = true;
noDataCount++;
}
else
{
if(currentValue > firstValue)
{
firstValue = currentValue;
highestPosition = currentPosition;
}
}
currentPosition++;
}
}
//the ArcGIS implementation uses 1 for first position value instead of 0 as in standard c++
return static_cast<int>( std::max_element( cellValueStack.begin(), cellValueStack.end() ) - cellValueStack.begin() ) + 1; //therefore... index + 1
return ++highestPosition; //therefore ++
}

///@endcond
Expand Down
8 changes: 4 additions & 4 deletions src/analysis/processing/qgsalgorithmrasterstackposition.h
Expand Up @@ -39,12 +39,12 @@ class QgsRasterStackPositionAlgorithmBase : public QgsProcessingAlgorithm
protected:
bool prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;
QVariantMap processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;
virtual int getPosition( std::vector<QgsRasterBlock>inputBlocks ) = 0;
virtual int findPosition( std::vector< std::unique_ptr< QgsRasterBlock > > &rasterBlockStack, int &row, int &col, bool &noDataInRasterBlockStack ) = 0;
double mNoDataValue = -9999;

private:
std::vector< QgsRasterAnalysisUtils::RasterLogicInput > mInputs;
bool mIgnoreNoData;
double mNoDataValue = -9999;
int mLayerWidth;
int mLayerHeight;
QgsRectangle mExtent;
Expand All @@ -64,7 +64,7 @@ class QgsRasterStackLowestPositionAlgorithm : public QgsRasterStackPositionAlgor
QgsRasterStackLowestPositionAlgorithm *createInstance() const override SIP_FACTORY;

protected:
int getPosition( std::vector<QgsRasterBlock>inputBlocks ) override;
int findPosition( std::vector< std::unique_ptr< QgsRasterBlock > > &rasterBlockStack, int &row, int &col, bool &noDataInRasterBlockStack ) override;
};

class QgsRasterStackHighestPositionAlgorithm : public QgsRasterStackPositionAlgorithmBase
Expand All @@ -78,7 +78,7 @@ class QgsRasterStackHighestPositionAlgorithm : public QgsRasterStackPositionAlgo
QgsRasterStackHighestPositionAlgorithm *createInstance() const override SIP_FACTORY;

protected:
int getPosition( std::vector<QgsRasterBlock>inputBlocks ) override;
int findPosition( std::vector< std::unique_ptr< QgsRasterBlock > > &rasterBlockStack, int &row, int &col, bool &noDataInRasterBlockStack ) override;
};

///@endcond PRIVATE
Expand Down

0 comments on commit 58629cb

Please sign in to comment.