QGIS API Documentation 3.39.0-Master (67e056379ed)
Loading...
Searching...
No Matches
qgsalgorithmzonalhistogram.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsalgorithmzonalhistogram.cpp
3 ---------------------
4 begin : May, 2018
5 copyright : (C) 2018 by Mathieu Pellerin
6 email : nirvn dot asia at gmail dot com
7 ***************************************************************************/
8
9/***************************************************************************
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 ***************************************************************************/
17
20#include "qgslogger.h"
21
23
24QString QgsZonalHistogramAlgorithm::name() const
25{
26 return QStringLiteral( "zonalhistogram" );
27}
28
29QString QgsZonalHistogramAlgorithm::displayName() const
30{
31 return QObject::tr( "Zonal histogram" );
32}
33
34QStringList QgsZonalHistogramAlgorithm::tags() const
35{
36 return QObject::tr( "raster,unique,values,count,area,statistics" ).split( ',' );
37}
38
39QString QgsZonalHistogramAlgorithm::group() const
40{
41 return QObject::tr( "Raster analysis" );
42}
43
44QString QgsZonalHistogramAlgorithm::groupId() const
45{
46 return QStringLiteral( "rasteranalysis" );
47}
48
49void QgsZonalHistogramAlgorithm::initAlgorithm( const QVariantMap & )
50{
51 addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT_RASTER" ),
52 QObject::tr( "Raster layer" ) ) );
53 addParameter( new QgsProcessingParameterBand( QStringLiteral( "RASTER_BAND" ),
54 QObject::tr( "Band number" ), 1, QStringLiteral( "INPUT_RASTER" ) ) );
55
56 addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT_VECTOR" ),
57 QObject::tr( "Vector layer containing zones" ), QList< int >() << static_cast< int >( Qgis::ProcessingSourceType::VectorPolygon ) ) );
58
59 addParameter( new QgsProcessingParameterString( QStringLiteral( "COLUMN_PREFIX" ), QObject::tr( "Output column prefix" ), QStringLiteral( "HISTO_" ), false, true ) );
60
61 addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Output zones" ), Qgis::ProcessingSourceType::VectorPolygon ) );
62}
63
64QString QgsZonalHistogramAlgorithm::shortHelpString() const
65{
66 return QObject::tr( "This algorithm appends fields representing counts of each unique value from a raster layer contained within zones defined as polygons." );
67}
68
69QgsZonalHistogramAlgorithm *QgsZonalHistogramAlgorithm::createInstance() const
70{
71 return new QgsZonalHistogramAlgorithm();
72}
73
74bool QgsZonalHistogramAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
75{
76 QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT_RASTER" ), context );
77 if ( !layer )
78 throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT_RASTER" ) ) );
79
80 mRasterBand = parameterAsInt( parameters, QStringLiteral( "RASTER_BAND" ), context );
81 mHasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( mRasterBand );
82 mNodataValue = layer->dataProvider()->sourceNoDataValue( mRasterBand );
83 mRasterInterface.reset( layer->dataProvider()->clone() );
84 mRasterExtent = layer->extent();
85 mCrs = layer->crs();
86 mCellSizeX = std::abs( layer->rasterUnitsPerPixelX() );
87 mCellSizeY = std::abs( layer->rasterUnitsPerPixelX() );
88 mNbCellsXProvider = mRasterInterface->xSize();
89 mNbCellsYProvider = mRasterInterface->ySize();
90
91 return true;
92}
93
94QVariantMap QgsZonalHistogramAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
95{
96
97 std::unique_ptr< QgsFeatureSource > zones( parameterAsSource( parameters, QStringLiteral( "INPUT_VECTOR" ), context ) );
98 if ( !zones )
99 throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT_VECTOR" ) ) );
100
101 const long count = zones->featureCount();
102 const double step = count > 0 ? 100.0 / count : 1;
103 long current = 0;
104
105 QList< double > uniqueValues;
106 QMap< QgsFeatureId, QHash< double, qgssize > > featuresUniqueValues;
107
108 // First loop through the zones to build up a list of unique values across all zones to determine sink fields list
109 QgsFeatureRequest request;
110 request.setNoAttributes();
111 if ( zones->sourceCrs() != mCrs )
112 {
113 request.setDestinationCrs( mCrs, context.transformContext() );
114 }
115 QgsFeatureIterator it = zones->getFeatures( request );
116 QgsFeature f;
117 while ( it.nextFeature( f ) )
118 {
119 if ( feedback->isCanceled() )
120 {
121 break;
122 }
123 feedback->setProgress( current * step );
124
125 if ( !f.hasGeometry() )
126 {
127 current++;
128 continue;
129 }
130
131 const QgsGeometry featureGeometry = f.geometry();
132 const QgsRectangle featureRect = featureGeometry.boundingBox().intersect( mRasterExtent );
133 if ( featureRect.isEmpty() )
134 {
135 current++;
136 continue;
137 }
138
139 int nCellsX, nCellsY;
140 QgsRectangle rasterBlockExtent;
141 QgsRasterAnalysisUtils::cellInfoForBBox( mRasterExtent, featureRect, mCellSizeX, mCellSizeY, nCellsX, nCellsY, mNbCellsXProvider, mNbCellsYProvider, rasterBlockExtent );
142
143 QHash< double, qgssize > fUniqueValues;
144 QgsRasterAnalysisUtils::statisticsFromMiddlePointTest( mRasterInterface.get(), mRasterBand, featureGeometry, nCellsX, nCellsY, mCellSizeX, mCellSizeY,
145 rasterBlockExtent, [ &fUniqueValues]( double value ) { fUniqueValues[value]++; }, false );
146
147 if ( fUniqueValues.count() < 1 )
148 {
149 // The cell resolution is probably larger than the polygon area. We switch to slower precise pixel - polygon intersection in this case
150 // TODO: eventually deal with weight if needed
151 QgsRasterAnalysisUtils::statisticsFromPreciseIntersection( mRasterInterface.get(), mRasterBand, featureGeometry, nCellsX, nCellsY, mCellSizeX, mCellSizeY,
152 rasterBlockExtent, [ &fUniqueValues]( double value, double ) { fUniqueValues[value]++; }, false );
153 }
154
155 for ( auto it = fUniqueValues.constBegin(); it != fUniqueValues.constEnd(); ++it )
156 {
157 if ( uniqueValues.indexOf( it.key() ) == -1 )
158 {
159 uniqueValues << it.key();
160 }
161 featuresUniqueValues[f.id()][it.key()] += it.value();
162 }
163
164 current++;
165 }
166
167 std::sort( uniqueValues.begin(), uniqueValues.end() );
168
169 const QString fieldPrefix = parameterAsString( parameters, QStringLiteral( "COLUMN_PREFIX" ), context );
170 QgsFields newFields;
171 for ( auto it = uniqueValues.constBegin(); it != uniqueValues.constEnd(); ++it )
172 {
173 newFields.append( QgsField( QStringLiteral( "%1%2" ).arg( fieldPrefix, mHasNoDataValue && *it == mNodataValue ? QStringLiteral( "NODATA" ) : QString::number( *it ) ), QMetaType::Type::LongLong, QString(), -1, 0 ) );
174 }
175 const QgsFields fields = QgsProcessingUtils::combineFields( zones->fields(), newFields );
176
177 QString dest;
178 std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, fields,
179 zones->wkbType(), zones->sourceCrs() ) );
180 if ( !sink )
181 throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );
182
183 it = zones->getFeatures( QgsFeatureRequest() );
184 while ( it.nextFeature( f ) )
185 {
186 QgsAttributes attributes = f.attributes();
187 const QHash< double, qgssize > fUniqueValues = featuresUniqueValues.value( f.id() );
188 for ( auto it = uniqueValues.constBegin(); it != uniqueValues.constEnd(); ++it )
189 {
190 attributes += fUniqueValues.value( *it, 0 );
191 }
192
193 QgsFeature outputFeature;
194 outputFeature.setGeometry( f.geometry() );
195 outputFeature.setAttributes( attributes );
196
197 if ( !sink->addFeature( outputFeature, QgsFeatureSink::FastInsert ) )
198 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
199 }
200
201 QVariantMap outputs;
202 outputs.insert( QStringLiteral( "OUTPUT" ), dest );
203 return outputs;
204}
205
207
208
209
@ VectorPolygon
Vector polygon layers.
A vector of attributes.
Wrapper for iterator of features from vector data provider or vector layer.
bool nextFeature(QgsFeature &f)
Fetch next feature and stores in f, returns true on success.
This class wraps a request for features to a vector layer (or directly its vector data provider).
QgsFeatureRequest & setDestinationCrs(const QgsCoordinateReferenceSystem &crs, const QgsCoordinateTransformContext &context)
Sets the destination crs for feature's geometries.
QgsFeatureRequest & setNoAttributes()
Set that no attributes will be fetched.
@ FastInsert
Use faster inserts, at the cost of updating the passed features to reflect changes made at the provid...
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition qgsfeature.h:58
QgsAttributes attributes
Definition qgsfeature.h:67
QgsFeatureId id
Definition qgsfeature.h:66
void setAttributes(const QgsAttributes &attrs)
Sets the feature's attributes.
QgsGeometry geometry
Definition qgsfeature.h:69
bool hasGeometry() const
Returns true if the feature has an associated geometry.
void setGeometry(const QgsGeometry &geometry)
Set the feature's geometry.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition qgsfeedback.h:53
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition qgsfeedback.h:61
Encapsulate a field in an attribute table or data source.
Definition qgsfield.h:53
Container of fields for a vector layer.
Definition qgsfields.h:46
bool append(const QgsField &field, Qgis::FieldOrigin origin=Qgis::FieldOrigin::Provider, int originIndex=-1)
Appends a field.
Definition qgsfields.cpp:69
A geometry is the spatial representation of a feature.
QgsRectangle boundingBox() const
Returns the bounding box of the geometry.
virtual QgsRectangle extent() const
Returns the extent of the layer.
QgsCoordinateReferenceSystem crs
Definition qgsmaplayer.h:82
Contains information about the context in which a processing algorithm is executed.
QgsCoordinateTransformContext transformContext() const
Returns the coordinate transform context.
Custom exception class for processing related exceptions.
Base class for providing feedback from a processing algorithm.
A raster band parameter for Processing algorithms.
A feature sink output for processing algorithms.
An input feature source (such as vector layers) parameter for processing algorithms.
A raster layer parameter for processing algorithms.
A string parameter for processing algorithms.
static QgsFields combineFields(const QgsFields &fieldsA, const QgsFields &fieldsB, const QString &fieldsBPrefix=QString())
Combines two field lists, avoiding duplicate field names (in a case-insensitive manner).
QgsRasterDataProvider * clone() const override=0
Clone itself, create deep copy.
virtual bool sourceHasNoDataValue(int bandNo) const
Returns true if source band has no data value.
virtual double sourceNoDataValue(int bandNo) const
Value representing no data value.
Represents a raster layer.
double rasterUnitsPerPixelX() const
Returns the number of raster units per each raster pixel in X axis.
QgsRasterDataProvider * dataProvider() override
Returns the source data provider.
A rectangle specified with double values.
bool isEmpty() const
Returns true if the rectangle has no area.
QgsRectangle intersect(const QgsRectangle &rect) const
Returns the intersection with the given rectangle.