Skip to content
Permalink
Browse files

add native Equal to frequency algorithm

  • Loading branch information
root676 authored and nyalldawson committed Jul 28, 2020
1 parent 21a7962 commit 175f749a11be8155ee3a2dacff3101094c05669a
@@ -127,6 +127,7 @@ SET(QGIS_ANALYSIS_SRCS
processing/qgsalgorithmrandompointsinpolygons.cpp
processing/qgsalgorithmrandompointsonlines.cpp
processing/qgsalgorithmrandomraster.cpp
processing/qgsalgorithmrasterequaltofrequency.cpp
processing/qgsalgorithmrasterlayeruniquevalues.cpp
processing/qgsalgorithmrasterlogicalop.cpp
processing/qgsalgorithmrasterize.cpp
@@ -0,0 +1,252 @@
/***************************************************************************
qgsalgorithmrasterequaltofrequency.cpp
---------------------
begin : June 2020
copyright : (C) 2020 by Clemens Raffler
email : clemens dot raffler at gmail dot com
***************************************************************************/

/***************************************************************************
* *
* 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 "qgsalgorithmrasterequaltofrequency.h"
#include "qgsrasterprojector.h"
#include "qgsrasterfilewriter.h"
#include "qgsrasteranalysisutils.h"

///@cond PRIVATE

QString QgsRasterEqualToFrequencyAlgorithm::displayName() const
{
return QObject::tr( "Equal to frequency" );
}

QString QgsRasterEqualToFrequencyAlgorithm::name() const
{
return QObject::tr( "equaltofrequency" );
}

QStringList QgsRasterEqualToFrequencyAlgorithm::tags() const
{
return QObject::tr( "cell,equal,frequency,pixel,stack" ).split( ',' );
}

QString QgsRasterEqualToFrequencyAlgorithm::group() const
{
return QObject::tr( "Raster analysis" );
}

QString QgsRasterEqualToFrequencyAlgorithm::groupId() const
{
return QStringLiteral( "rasteranalysis" );
}

QString QgsRasterEqualToFrequencyAlgorithm::shortHelpString() const
{
return QObject::tr( "The Equal to frequency algorithm evaluates on a cell-by-cell basis the frequency "
"(number of times) the values of an input stack of rasters are equal "
"to the value of an input raster. \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 - use GDAL to use other bands in the analysis."
"The input value layer serves as reference layer for the sample layers."
"Any NoData cells in the value raster or the data layer stack will result in a NoData cell "
"in the output raster if the ignore NoData parameter is not checked."
"The output NoData value can be set manually. The output rasters extent and resolution "
"is defined by the input raster layer and is always of int32 type." );
}

QgsRasterEqualToFrequencyAlgorithm *QgsRasterEqualToFrequencyAlgorithm::createInstance() const
{
return new QgsRasterEqualToFrequencyAlgorithm();
}

void QgsRasterEqualToFrequencyAlgorithm::initAlgorithm( const QVariantMap & )
{
addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral("INPUT_VALUE_RASTER"), QObject::tr("Input value raster") ) );
addParameter( new QgsProcessingParameterBand( QStringLiteral("INPUT_VALUE_RASTER_BAND"), QObject::tr("Value raster band"), 1, QStringLiteral("INPUT_VALUE_RASTER") ) );


addParameter( new QgsProcessingParameterMultipleLayers( QStringLiteral( "INPUT_RASTERS" ),
QObject::tr( "Input raster layers" ), QgsProcessing::TypeRaster ) );

addParameter( new QgsProcessingParameterBoolean( QStringLiteral( "IGNORE_NODATA" ), QObject::tr( "Ignore NoData values" ), false ) );

std::unique_ptr< QgsProcessingParameterNumber > output_nodata_parameter = qgis::make_unique< QgsProcessingParameterNumber >( QStringLiteral( "OUTPUT_NODATA_VALUE" ), QObject::tr( "Output NoData value" ), QgsProcessingParameterNumber::Double, -9999, true );
output_nodata_parameter->setFlags( output_nodata_parameter->flags() | QgsProcessingParameterDefinition::FlagAdvanced );
addParameter( output_nodata_parameter.release() );

addParameter( new QgsProcessingParameterRasterDestination( QStringLiteral( "OUTPUT" ),
QObject::tr( "Output layer" ) ) );
addOutput( new QgsProcessingOutputNumber( QStringLiteral( "EQUAL_VALUES_COUNT"), QObject::tr("Count of equal value occurrances") ) );
addOutput( new QgsProcessingOutputNumber( QStringLiteral( "FOUND_LOCATIONS_COUNT" ), QObject::tr("Count of cells with equal value occurrances") ) );
addOutput( new QgsProcessingOutputNumber( QStringLiteral( "MEAN_FREQUENCY_PER_LOCATION" ), QObject::tr("Mean frequency at valid cell locations" ) ) );
addOutput( new QgsProcessingOutputString( QStringLiteral( "EXTENT" ), QObject::tr( "Extent" ) ) );
addOutput( new QgsProcessingOutputString( QStringLiteral( "CRS_AUTHID" ), QObject::tr( "CRS authority identifier" ) ) );
addOutput( new QgsProcessingOutputNumber( QStringLiteral( "WIDTH_IN_PIXELS" ), QObject::tr( "Width in pixels" ) ) );
addOutput( new QgsProcessingOutputNumber( QStringLiteral( "HEIGHT_IN_PIXELS" ), QObject::tr( "Height in pixels" ) ) );
addOutput( new QgsProcessingOutputNumber( QStringLiteral( "TOTAL_PIXEL_COUNT" ), QObject::tr( "Total pixel count" ) ) );
}

bool QgsRasterEqualToFrequencyAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
{
QgsRasterLayer *inputValueRaster = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT_VALUE_RASTER" ), context );
if ( !inputValueRaster )
throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT_VALUE_RASTER" ) ) );

mInputValueRasterBand = parameterAsInt( parameters, QStringLiteral( "INPUT_VALUE_RASTER_BAND"), context );
mIgnoreNoData = parameterAsBool( parameters, QStringLiteral( "IGNORE_NODATA" ), context );

mInputValueRasterInterface.reset( inputValueRaster->dataProvider()->clone() );
mNoDataValue = parameterAsDouble( parameters, QStringLiteral( "OUTPUT_NODATA_VALUE" ), context );
mCrs = inputValueRaster->crs();
mRasterUnitsPerPixelX = inputValueRaster->rasterUnitsPerPixelX();
mRasterUnitsPerPixelY = inputValueRaster->rasterUnitsPerPixelY();
mLayerWidth = inputValueRaster->width();
mLayerHeight = inputValueRaster->height();
mExtent = inputValueRaster->extent();

const QList< QgsMapLayer * > layers = parameterAsLayerList( parameters, QStringLiteral( "INPUT_RASTERS" ), context );
QList< QgsRasterLayer * > rasterLayers;
rasterLayers.reserve( layers.count() );
for ( QgsMapLayer *l : layers )
{
if ( feedback->isCanceled() )
break; //in case some slow data sources are loaded

if ( l->type() == QgsMapLayerType::RasterLayer )
{
QgsRasterLayer *layer = qobject_cast< QgsRasterLayer * >( l );
QgsRasterAnalysisUtils::RasterLogicInput input;
const int band = 1; //could be made dynamic
input.hasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( band );
input.sourceDataProvider.reset( layer->dataProvider()->clone() );
input.interface = input.sourceDataProvider.get();
// add projector if necessary
if ( layer->crs() != mCrs )
{
input.projector = qgis::make_unique< QgsRasterProjector >();
input.projector->setInput( input.sourceDataProvider.get() );
input.projector->setCrs( layer->crs(), mCrs, context.transformContext() );
input.interface = input.projector.get();
}
mInputs.emplace_back( std::move( input ) );
}
}

return true;
}

QVariantMap QgsRasterEqualToFrequencyAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
{
const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context );
QFileInfo fi( outputFile );
const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );

std::unique_ptr< QgsRasterFileWriter > writer = qgis::make_unique< QgsRasterFileWriter >( outputFile );
writer->setOutputProviderKey( QStringLiteral( "gdal" ) );
writer->setOutputFormat( outputFormat );
std::unique_ptr<QgsRasterDataProvider > provider( writer->createOneBandRaster( Qgis::Int32, mLayerWidth, mLayerHeight, mExtent, mCrs ) );
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, mNoDataValue );
qgssize layerSize = static_cast< qgssize >( mLayerWidth ) * static_cast< qgssize >( mLayerHeight );

int maxWidth = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_WIDTH;
int maxHeight = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_HEIGHT;
int nbBlocksWidth = static_cast< int>( std::ceil( 1.0 * mLayerWidth / maxWidth ) );
int nbBlocksHeight = static_cast< int >( std::ceil( 1.0 * mLayerHeight / maxHeight ) );
int nbBlocks = nbBlocksWidth * nbBlocksHeight;
provider->setEditable( true );

QgsRasterIterator iter( mInputValueRasterInterface.get() );
iter.startRasterRead( mInputValueRasterBand, mLayerWidth, mLayerHeight, mExtent );
int iterLeft = 0;
int iterTop = 0;
int iterCols = 0;
int iterRows = 0;
QgsRectangle blockExtent;

unsigned long long equalValuesCount = 0;
unsigned long long noDataLocationsCount = 0;
std::unique_ptr< QgsRasterBlock > inputBlock;
while ( iter.readNextRasterPart( 1, iterCols, iterRows, inputBlock, iterLeft, iterTop, &blockExtent ) )
{
std::vector< std::unique_ptr< QgsRasterBlock > > inputBlocks;
for ( const QgsRasterAnalysisUtils::RasterLogicInput &i : mInputs )
{
if ( feedback->isCanceled() )
break; //in case some slow data sources are loaded
for ( int band : i.bands )
{
if ( feedback->isCanceled() )
break; //in case some slow data sources are loaded
std::unique_ptr< QgsRasterBlock > b( i.interface->block( band, blockExtent, iterCols, iterRows ) );
inputBlocks.emplace_back( std::move( b ) );
}
}

std::unique_ptr< QgsRasterBlock > outputBlock = qgis::make_unique<QgsRasterBlock>( Qgis::Int32, iterCols, iterRows);
feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
for ( int row = 0; row < iterRows; row++ )
{
if ( feedback->isCanceled() )
break;

for ( int col = 0; col < iterCols; col++ )
{
bool noDataInStack = false;
std::vector<double> cellValues = QgsRasterAnalysisUtils::getCellValuesFromBlockStack( inputBlocks, row, col, noDataInStack );

bool valueRasterCellIsNoData = false;
double searchValue = inputBlock->valueAndNoData( row, col, valueRasterCellIsNoData );

if ( (valueRasterCellIsNoData || noDataInStack) && !mIgnoreNoData )
{
//output cell will always be NoData if NoData occurs in valueRaster or cellValueStack and NoData is not ignored
//this saves unnecessary iterations on the cellValueStack
outputBlock->setValue( row, col, mNoDataValue );
noDataLocationsCount++;
}
else
{
int equalCount = static_cast<int>(std::count(cellValues.begin(), cellValues.end(), searchValue));
outputBlock->setValue( row, col, equalCount );
equalValuesCount += equalCount;
}
}
}
provider->writeBlock( outputBlock.get(), 1, iterLeft, iterTop );
}
provider->setEditable( false );

unsigned long long foundLocationsCount = layerSize - noDataLocationsCount;
double meanEqualCountPerValidLocation = static_cast<double>(equalValuesCount) / static_cast<double>(foundLocationsCount * mInputs.size() );

QVariantMap outputs;
outputs.insert( QStringLiteral( "EQUAL_VALUES_COUNT" ), equalValuesCount);
outputs.insert( QStringLiteral( "FOUND_LOCATIONS_COUNT" ), foundLocationsCount);
outputs.insert( QStringLiteral( "MEAN_FREQUENCY_PER_LOCATION" ), meanEqualCountPerValidLocation );
outputs.insert( QStringLiteral( "EXTENT" ), mExtent.toString() );
outputs.insert( QStringLiteral( "CRS_AUTHID" ), mCrs.authid() );
outputs.insert( QStringLiteral( "WIDTH_IN_PIXELS" ), mLayerWidth );
outputs.insert( QStringLiteral( "HEIGHT_IN_PIXELS" ), mLayerHeight );
outputs.insert( QStringLiteral( "TOTAL_PIXEL_COUNT" ), layerSize );
outputs.insert( QStringLiteral( "OUTPUT" ), outputFile );

return outputs;
}



///@endcond


@@ -0,0 +1,68 @@
/***************************************************************************
qgsalgorithmrasterequaltofrequency.h
---------------------
begin : June 2020
copyright : (C) 2020 by Clemens Raffler
email : clemens dot raffler at gmail dot com
***************************************************************************/

/***************************************************************************
* *
* 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 QGSALGORITHMRASTEREQUALTOFREQUENCY_H
#define QGSALGORITHMRASTEREQUALTOFREQUENCY_H

#define SIP_NO_FILE

#include "qgis_sip.h"
#include "qgsapplication.h"
#include "qgsprocessingalgorithm.h"
#include "qgsrasterprojector.h"
#include "qgsrasteranalysisutils.h"

///@cond PRIVATE

class QgsRasterEqualToFrequencyAlgorithm : public QgsProcessingAlgorithm
{

public:
QgsRasterEqualToFrequencyAlgorithm() = 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;
QgsRasterEqualToFrequencyAlgorithm *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;

private:
std::unique_ptr< QgsRasterInterface > mInputValueRasterInterface;
int mInputValueRasterBand;
std::vector< QgsRasterAnalysisUtils::RasterLogicInput > mInputs;
bool mIgnoreNoData;
double mNoDataValue = -9999;
int mLayerWidth;
int mLayerHeight;
QgsRectangle mExtent;
QgsCoordinateReferenceSystem mCrs;
double mRasterUnitsPerPixelX;
double mRasterUnitsPerPixelY;
};

///@endcond PRIVATE

#endif // QGSALGORITHMRASTEREQUALTOFREQUENCY_H

@@ -122,6 +122,7 @@
#include "qgsalgorithmrandompointsinpolygons.h"
#include "qgsalgorithmrandompointsonlines.h"
#include "qgsalgorithmrandomraster.h"
#include "qgsalgorithmrasterequaltofrequency.h"
#include "qgsalgorithmrasterlayeruniquevalues.h"
#include "qgsalgorithmrasterlogicalop.h"
#include "qgsalgorithmrasterize.h"
@@ -355,6 +356,7 @@ void QgsNativeAlgorithms::loadAlgorithms()
addAlgorithm( new QgsRandomPointsOnLinesAlgorithm() );
addAlgorithm( new QgsRandomPoissonRasterAlgorithm() );
addAlgorithm( new QgsRandomUniformRasterAlgorithm() );
addAlgorithm( new QgsRasterEqualToFrequencyAlgorithm() );
addAlgorithm( new QgsRasterLayerUniqueValuesReportAlgorithm() );
addAlgorithm( new QgsRasterLayerZonalStatsAlgorithm() );
addAlgorithm( new QgsRasterLogicalAndAlgorithm() );

0 comments on commit 175f749

Please sign in to comment.
You can’t perform that action at this time.