Skip to content

Commit

Permalink
Add new algorithm to calculate zonal statistics into a new sink
Browse files Browse the repository at this point in the history
  • Loading branch information
m-kuhn committed Sep 9, 2020
1 parent c20ebcf commit 94fe2f8
Show file tree
Hide file tree
Showing 8 changed files with 447 additions and 267 deletions.
27 changes: 10 additions & 17 deletions python/analysis/auto_generated/vector/qgszonalstatistics.sip.in
Expand Up @@ -107,23 +107,6 @@ added to ``polygonLayer`` for each statistic calculated.
%End


QgsZonalStatistics( QgsFeatureSource *source,
QgsFeatureSink *sink,
QgsRasterInterface *rasterInterface,
const QgsCoordinateReferenceSystem &rasterCrs,
const QMap<QgsZonalStatistics::Statistic, int> &statFieldIndexes,
const QgsFields &fields,
double rasterUnitsPerPixelX,
double rasterUnitsPerPixelY,
int rasterBand = 1,
QgsZonalStatistics::Statistics stats = QgsZonalStatistics::Statistic::All );
%Docstring


.. versionadded:: 3.16
%End


QgsZonalStatistics::Result calculateStatistics( QgsFeedback *feedback );
%Docstring
Runs the calculation.
Expand All @@ -145,6 +128,16 @@ Returns a short, friendly display name for a ``statistic``, suitable for use in
.. seealso:: :py:func:`displayName`

.. versionadded:: 3.12
%End

static QMap<QgsZonalStatistics::Statistic, QVariant> calculateStatistics( QgsRasterInterface *rasterInterface, const QgsGeometry &geometry, int cellSizeX, int cellSizeY, int rasterBand, QgsZonalStatistics::Statistics statistics );
%Docstring
Calculates the specified ``statistics`` for the pixels of ``rasterBand``
in ``rasterInterface`` (a raster layer :py:func:`~QgsZonalStatistics.dataProvider` ) within polygon ``geometry``.

Returns a map of statistic to result value.

.. versionadded:: 3.16
%End

public:
Expand Down
1 change: 1 addition & 0 deletions src/analysis/CMakeLists.txt
Expand Up @@ -197,6 +197,7 @@ SET(QGIS_ANALYSIS_SRCS
processing/qgsalgorithmwritevectortiles.cpp
processing/qgsalgorithmzonalhistogram.cpp
processing/qgsalgorithmzonalstatistics.cpp
processing/qgsalgorithmzonalstatisticsfeaturebased.cpp
processing/qgsbookmarkalgorithms.cpp
processing/qgsprojectstylealgorithms.cpp
processing/qgsstylealgorithms.cpp
Expand Down
56 changes: 9 additions & 47 deletions src/analysis/processing/qgsalgorithmzonalstatistics.cpp
Expand Up @@ -64,7 +64,7 @@ QString QgsZonalStatisticsAlgorithm::groupId() const
QString QgsZonalStatisticsAlgorithm::shortHelpString() const
{
return QObject::tr( "This algorithm calculates statistics of a raster layer for each feature "
"of an overlapping polygon vector layer." );
"of an overlapping polygon vector layer. The reults will be written in place." );
}

QgsProcessingAlgorithm::Flags QgsZonalStatisticsAlgorithm::flags() const
Expand All @@ -89,14 +89,14 @@ void QgsZonalStatisticsAlgorithm::initAlgorithm( const QVariantMap & )
addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT_RASTER" ), QObject::tr( "Raster layer" ) ) );
addParameter( new QgsProcessingParameterBand( QStringLiteral( "RASTER_BAND" ),
QObject::tr( "Raster band" ), 1, QStringLiteral( "INPUT_RASTER" ) ) );
addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT_VECTOR" ), QObject::tr( "Feature source containing zones" ),
addParameter( new QgsProcessingParameterVectorLayer( QStringLiteral( "INPUT_VECTOR" ), QObject::tr( "Vector layer containing zones" ),
QList< int >() << QgsProcessing::TypeVectorPolygon ) );
addParameter( new QgsProcessingParameterString( QStringLiteral( "COLUMN_PREFIX" ), QObject::tr( "Output column prefix" ), QStringLiteral( "_" ) ) );

addParameter( new QgsProcessingParameterEnum( QStringLiteral( "STATISTICS" ), QObject::tr( "Statistics to calculate" ),
statChoices, true, QVariantList() << 0 << 1 << 2 ) );

addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Zonal statistics" ), QgsProcessing::TypeVectorPolygon, QVariant(), false, true, true ) );
addOutput( new QgsProcessingOutputVectorLayer( QStringLiteral( "INPUT_VECTOR" ), QObject::tr( "Zonal statistics" ), QgsProcessing::TypeVectorPolygon ) );
}

bool QgsZonalStatisticsAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
Expand Down Expand Up @@ -129,62 +129,24 @@ bool QgsZonalStatisticsAlgorithm::prepareAlgorithm( const QVariantMap &parameter

QVariantMap QgsZonalStatisticsAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
{
std::unique_ptr<QgsFeatureSource> source( parameterAsSource( parameters, QStringLiteral( "INPUT_VECTOR" ), context ) );
if ( !source )
QgsVectorLayer *layer = parameterAsVectorLayer( parameters, QStringLiteral( "INPUT_VECTOR" ), context );
if ( !layer )
throw QgsProcessingException( QObject::tr( "Invalid zones layer" ) );

QgsFields fields = source->fields();
QMap<QgsZonalStatistics::Statistic, int> statFieldIndexes;

for ( QgsZonalStatistics::Statistic stat :
{
QgsZonalStatistics::Count,
QgsZonalStatistics::Sum,
QgsZonalStatistics::Mean,
QgsZonalStatistics::Median,
QgsZonalStatistics::StDev,
QgsZonalStatistics::Min,
QgsZonalStatistics::Max,
QgsZonalStatistics::Range,
QgsZonalStatistics::Minority,
QgsZonalStatistics::Majority,
QgsZonalStatistics::Variety,
QgsZonalStatistics::Variance
} )
{
if ( mStats & stat )
{
QgsField field = QgsField( mPrefix + QgsZonalStatistics::shortName( stat ), QVariant::Double, QStringLiteral( "double precision" ) );
if ( fields.names().contains( field.name() ) )
{
throw QgsProcessingException( QObject::tr( "Field %1 already exists" ).arg( field.name() ) );
}
fields.append( field );
statFieldIndexes.insert( stat, fields.count() - 1 );
}
}

QString dest;
std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, fields, QgsWkbTypes::MultiPolygon, source->sourceCrs() ) );
if ( !sink )
throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );

QgsZonalStatistics zs( source.get(),
sink.get(),
QgsZonalStatistics zs( layer,
mInterface.get(),
mCrs,
statFieldIndexes,
fields,
mPixelSizeX,
mPixelSizeY,
mPrefix,
mBand,
mStats
QgsZonalStatistics::Statistics( mStats )
);

zs.calculateStatistics( feedback );

QVariantMap outputs;
outputs.insert( QStringLiteral( "OUTPUT" ), dest );
outputs.insert( QStringLiteral( "INPUT_VECTOR" ), layer->id() );
return outputs;
}

Expand Down
189 changes: 189 additions & 0 deletions src/analysis/processing/qgsalgorithmzonalstatisticsfeaturebased.cpp
@@ -0,0 +1,189 @@
/***************************************************************************
qgsalgorithmzonalstatisticsfeaturebased.cpp
---------------------
begin : September 2020
copyright : (C) 2020 by Matthias Kuhn
email : matthias@opengis.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 "qgsalgorithmzonalstatisticsfeaturebased.h"

///@cond PRIVATE

const std::vector< QgsZonalStatistics::Statistic > STATS
{
QgsZonalStatistics::Count,
QgsZonalStatistics::Sum,
QgsZonalStatistics::Mean,
QgsZonalStatistics::Median,
QgsZonalStatistics::StDev,
QgsZonalStatistics::Min,
QgsZonalStatistics::Max,
QgsZonalStatistics::Range,
QgsZonalStatistics::Minority,
QgsZonalStatistics::Majority,
QgsZonalStatistics::Variety,
QgsZonalStatistics::Variance,
};

QString QgsZonalStatisticsFeatureBasedAlgorithm::name() const
{
return QStringLiteral( "zonalstatisticsfb" );
}

QString QgsZonalStatisticsFeatureBasedAlgorithm::displayName() const
{
return QObject::tr( "Zonal statistics (feature based)" );
}

QStringList QgsZonalStatisticsFeatureBasedAlgorithm::tags() const
{
return QObject::tr( "stats,statistics,zones,layer,sum,maximum,minimum,mean,count,standard,deviation,"
"median,range,majority,minority,variety,variance,summary,raster" ).split( ',' );
}

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

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

QString QgsZonalStatisticsFeatureBasedAlgorithm::shortHelpString() const
{
return QObject::tr( "This algorithm calculates statistics of a raster layer for each feature "
"of an overlapping polygon vector layer." );
}

QList<int> QgsZonalStatisticsFeatureBasedAlgorithm::inputLayerTypes() const
{
return QList<int>() << QgsProcessing::TypeVectorPolygon;
}

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

void QgsZonalStatisticsFeatureBasedAlgorithm::initParameters( const QVariantMap &configuration )
{
Q_UNUSED( configuration )
QStringList statChoices;
statChoices.reserve( STATS.size() );
for ( QgsZonalStatistics::Statistic stat : STATS )
{
statChoices << QgsZonalStatistics::displayName( stat );
}

addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT_RASTER" ), QObject::tr( "Raster layer" ) ) );
addParameter( new QgsProcessingParameterBand( QStringLiteral( "RASTER_BAND" ),
QObject::tr( "Raster band" ), 1, QStringLiteral( "INPUT_RASTER" ) ) );

addParameter( new QgsProcessingParameterString( QStringLiteral( "COLUMN_PREFIX" ), QObject::tr( "Output column prefix" ), QStringLiteral( "_" ) ) );

addParameter( new QgsProcessingParameterEnum( QStringLiteral( "STATISTICS" ), QObject::tr( "Statistics to calculate" ),
statChoices, true, QVariantList() << 0 << 1 << 2 ) );
}

QString QgsZonalStatisticsFeatureBasedAlgorithm::outputName() const
{
return QObject::tr( "Zonal Statistics" );
}

QgsFields QgsZonalStatisticsFeatureBasedAlgorithm::outputFields( const QgsFields &inputFields ) const
{
Q_UNUSED( inputFields )
return mOutputFields;
}

QgsProcessingFeatureSource::Flag QgsZonalStatisticsFeatureBasedAlgorithm::sourceFlags() const
{
return QgsProcessingFeatureSource::FlagSkipGeometryValidityChecks;
}

bool QgsZonalStatisticsFeatureBasedAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
{
mPrefix = parameterAsString( parameters, QStringLiteral( "COLUMN_PREFIX" ), context );

const QList< int > stats = parameterAsEnums( parameters, QStringLiteral( "STATISTICS" ), context );
mStats = nullptr;
for ( int s : stats )
{
mStats |= STATS.at( s );
}

QgsRasterLayer *rasterLayer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT_RASTER" ), context );
if ( !rasterLayer )
throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT_RASTER" ) ) );

mBand = parameterAsInt( parameters, QStringLiteral( "RASTER_BAND" ), context );
if ( mBand < 1 || mBand > rasterLayer->bandCount() )
throw QgsProcessingException( QObject::tr( "Invalid band number for BAND (%1): Valid values for input raster are 1 to %2" ).arg( mBand )
.arg( rasterLayer->bandCount() ) );

if ( !rasterLayer->dataProvider() )
throw QgsProcessingException( QObject::tr( "Invalid raster layer. Layer %1 is invalid." ).arg( rasterLayer->id() ) );

mRaster.reset( rasterLayer->dataProvider()->clone() );
mCrs = rasterLayer->crs();
mPixelSizeX = rasterLayer->rasterUnitsPerPixelX();
mPixelSizeY = rasterLayer->rasterUnitsPerPixelY();
QgsFeatureSource *source = parameterAsSource( parameters, inputParameterName(), context );

mOutputFields = source->fields();

for ( QgsZonalStatistics::Statistic stat : STATS )
{
if ( mStats & stat )
{
QgsField field = QgsField( mPrefix + QgsZonalStatistics::shortName( stat ), QVariant::Double, QStringLiteral( "double precision" ) );
if ( mOutputFields.names().contains( field.name() ) )
{
throw QgsProcessingException( QObject::tr( "Field %1 already exists" ).arg( field.name() ) );
}
mOutputFields.append( field );
mStatFieldsMapping.insert( stat, mOutputFields.size() - 1 );
}
}

return true;
}

QgsFeatureList QgsZonalStatisticsFeatureBasedAlgorithm::processFeature( const QgsFeature &feature, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
{
Q_UNUSED( context )
Q_UNUSED( feedback )
QgsAttributes attributes = feature.attributes();
attributes.resize( mOutputFields.size() );

QMap<QgsZonalStatistics::Statistic, QVariant> results = QgsZonalStatistics::calculateStatistics( mRaster.get(), feature.geometry(), mPixelSizeX, mPixelSizeY, mBand, mStats );
for ( const auto &result : results.toStdMap() )
{
attributes.replace( mStatFieldsMapping.value( result.first ), result.second );
}

QgsFeature resultFeature = feature;
resultFeature.setAttributes( attributes );

return QgsFeatureList { resultFeature };
}

bool QgsZonalStatisticsFeatureBasedAlgorithm::supportInPlaceEdit( const QgsMapLayer *layer ) const
{
Q_UNUSED( layer )
return false;
}

///@endcond

0 comments on commit 94fe2f8

Please sign in to comment.