Skip to content
Permalink
Browse files

Merge pull request #38650 from m-kuhn/zonal_stats_new_layer

[processing] Add Zonal statistics algorithm that creates a new output
  • Loading branch information
m-kuhn committed Sep 10, 2020
2 parents 8639b1d + 8f7abb5 commit 8b7f600c33ae230bfc6ed6ae3e92c43516a8e1ab
@@ -42,6 +42,17 @@ A class that calculates raster statistics (count, sum, mean) for a polygon or mu
typedef QFlags<QgsZonalStatistics::Statistic> Statistics;


enum Result
{
Success,
LayerTypeWrong,
LayerInvalid,
RasterInvalid,
RasterBandInvalid,
FailedToCreateField,
Canceled
};

QgsZonalStatistics( QgsVectorLayer *polygonLayer,
QgsRasterLayer *rasterLayer,
const QString &attributePrefix = QString(),
@@ -95,11 +106,10 @@ added to ``polygonLayer`` for each statistic calculated.
.. versionadded:: 3.2
%End

int calculateStatistics( QgsFeedback *feedback );
%Docstring
Starts the calculation

:return: 0 in case of success
QgsZonalStatistics::Result calculateStatistics( QgsFeedback *feedback );
%Docstring
Runs the calculation.
%End

static QString displayName( QgsZonalStatistics::Statistic statistic );
@@ -118,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, double cellSizeX, double 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:
@@ -0,0 +1,9 @@
{
"type": "FeatureCollection",
"name": "stats",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [
{ "type": "Feature", "properties": { "fid": "polygon_mask.0", "stats_count": 24882.0, "stats_sum": 89005.0, "stats_mean": 3.5770838357045251, "stats_median": 6.0, "stats_stdev": 2.4739349627559846, "stats_min": 1.0, "stats_max": 6.0, "stats_range": 5.0, "stats_minority": 4.0, "stats_majority": 6.0, "stats_variety": 4.0, "stats_variance": 6.1203541999464539 }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 18.6744202640594, 45.798496472749385 ], [ 18.683204095342351, 45.807861586984885 ], [ 18.68966279481511, 45.80456765025378 ], [ 18.694829754393322, 45.800175734612303 ], [ 18.694442232424951, 45.796235927933921 ], [ 18.683849965289628, 45.790358511413707 ], [ 18.675970351932861, 45.790229337424257 ], [ 18.673193111159573, 45.789970989445344 ], [ 18.6744202640594, 45.798496472749385 ] ] ] ] } },
{ "type": "Feature", "properties": { "fid": "polygon_mask.1", "stats_count": 14055.0, "stats_sum": 42165.0, "stats_mean": 3.0, "stats_median": 3.0, "stats_stdev": 0.0, "stats_min": 3.0, "stats_max": 3.0, "stats_range": 0.0, "stats_minority": 3.0, "stats_majority": 3.0, "stats_variety": 1.0, "stats_variance": 0.0 }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 18.690696186730751, 45.786612465719514 ], [ 18.694442232424951, 45.791133555350441 ], [ 18.70090093189771, 45.785643660798598 ], [ 18.699738365992616, 45.781122571167664 ], [ 18.688435641915287, 45.777957808426017 ], [ 18.680491441563792, 45.779895418267841 ], [ 18.6800393326007, 45.783318528988403 ], [ 18.690696186730751, 45.786612465719514 ] ] ] ] } }
]
}
@@ -2163,4 +2163,33 @@ tests:
- '<SrcDataSource>.*\/multilines\.gml</SrcDataSource>'
- '<SrcLayer>multilines</SrcLayer>'

- algorithm: native:zonalstatisticsfb
name: Test (native:zonalstatisticsfb)
params:
COLUMN_PREFIX: stats_
INPUT:
name: custom/zonal_stats.shp
type: vector
INPUT_RASTER:
name: custom/dem_zones.tif
type: raster
RASTER_BAND: 1
STATISTICS:
- 0
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
results:
OUTPUT:
name: expected/stats.geojson
type: vector

# See ../README.md for a description of the file format
@@ -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
@@ -64,12 +64,12 @@ 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 results will be written in place." );
}

QgsProcessingAlgorithm::Flags QgsZonalStatisticsAlgorithm::flags() const
{
return QgsProcessingAlgorithm::flags() | QgsProcessingAlgorithm::FlagNoThreading;
return QgsProcessingAlgorithm::flags() | QgsProcessingAlgorithm::FlagNoThreading | QgsProcessingAlgorithm::FlagDeprecated;
}

QgsZonalStatisticsAlgorithm *QgsZonalStatisticsAlgorithm::createInstance() const
@@ -0,0 +1,184 @@
/***************************************************************************
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" );
}

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;
}

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();
std::unique_ptr<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 ( auto result = results.constBegin(); result != results.constEnd(); ++result )
{
attributes.replace( mStatFieldsMapping.value( result.key() ), result.value() );
}

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

return QgsFeatureList { resultFeature };
}

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

///@endcond
@@ -0,0 +1,73 @@
/***************************************************************************
qgsalgorithmzonalstatisticsfeaturebased.h
------------------------------
begin : September 2020
copyright : (C) 2020 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. *
* *
***************************************************************************/

#ifndef QGSALGORITHMZONALSTATISTICSFEATUREBASED_H
#define QGSALGORITHMZONALSTATISTICSFEATUREBASED_H

#define SIP_NO_FILE

#include "qgis_sip.h"
#include "qgsprocessingalgorithm.h"
#include "qgsvectorlayer.h"
#include "vector/qgszonalstatistics.h"

///@cond PRIVATE

/**
* Native zonal statistics algorithm.
*/
class QgsZonalStatisticsFeatureBasedAlgorithm : public QgsProcessingFeatureBasedAlgorithm
{

public:

QgsZonalStatisticsFeatureBasedAlgorithm() = default;
QString name() const override;
QString displayName() const override;
QStringList tags() const override;
QString group() const override;
QString groupId() const override;
QString shortHelpString() const override;
QList<int> inputLayerTypes() const override;

QgsZonalStatisticsFeatureBasedAlgorithm *createInstance() const override SIP_FACTORY;

protected:

void initParameters( const QVariantMap &configuration = QVariantMap() ) override;
QString outputName() const override;
QgsFields outputFields( const QgsFields &inputFields ) const override;

bool prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;
QgsFeatureList processFeature( const QgsFeature &feature, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;
bool supportInPlaceEdit( const QgsMapLayer *layer ) const override;

private:
std::unique_ptr< QgsRasterInterface > mRaster;
int mBand;
QString mPrefix;
QgsZonalStatistics::Statistics mStats = QgsZonalStatistics::All;
QgsCoordinateReferenceSystem mCrs;
double mPixelSizeX;
double mPixelSizeY;
QgsFields mOutputFields;
QMap<QgsZonalStatistics::Statistic, int> mStatFieldsMapping;
};

///@endcond PRIVATE

#endif // QGSALGORITHMZONALSTATISTICSFEATUREBASED_H
@@ -192,6 +192,7 @@
#include "qgsalgorithmwritevectortiles.h"
#include "qgsalgorithmzonalhistogram.h"
#include "qgsalgorithmzonalstatistics.h"
#include "qgsalgorithmzonalstatisticsfeaturebased.h"
#include "qgsalgorithmpolygonstolines.h"
#include "qgsbookmarkalgorithms.h"
#include "qgsprojectstylealgorithms.h"
@@ -441,6 +442,7 @@ void QgsNativeAlgorithms::loadAlgorithms()
addAlgorithm( new QgsWriteVectorTilesMbtilesAlgorithm() );
addAlgorithm( new QgsZonalHistogramAlgorithm() );
addAlgorithm( new QgsZonalStatisticsAlgorithm() );
addAlgorithm( new QgsZonalStatisticsFeatureBasedAlgorithm() );
addAlgorithm( new QgsPolygonsToLinesAlgorithm() );
addAlgorithm( new QgsDensifyGeometriesByIntervalAlgorithm() );
addAlgorithm( new QgsDensifyGeometriesByCountAlgorithm() );

0 comments on commit 8b7f600

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