QGIS API Documentation 3.41.0-Master (3c143d501a8)
Loading...
Searching...
No Matches
qgszonalstatistics.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgszonalstatistics.cpp - description
3 ----------------------------
4 begin : August 29th, 2009
5 copyright : (C) 2009 by Marco Hugentobler
6 email : marco at hugis dot net
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
18#include "qgszonalstatistics.h"
19
20#include "qgsfeatureiterator.h"
21#include "qgsfeedback.h"
22#include "qgsgeometry.h"
24#include "qgsvectorlayer.h"
26#include "qgsrasterlayer.h"
27#include "qgsproject.h"
28
29#include <QFile>
30
31QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer *polygonLayer, QgsRasterLayer *rasterLayer, const QString &attributePrefix, int rasterBand, Qgis::ZonalStatistics stats )
32 : QgsZonalStatistics( polygonLayer,
33 rasterLayer ? rasterLayer->dataProvider() : nullptr,
34 rasterLayer ? rasterLayer->crs() : QgsCoordinateReferenceSystem(),
35 rasterLayer ? rasterLayer->rasterUnitsPerPixelX() : 0,
36 rasterLayer ? rasterLayer->rasterUnitsPerPixelY() : 0,
37 attributePrefix,
38 rasterBand,
39 stats )
40{
41}
42
44 const QgsCoordinateReferenceSystem &rasterCrs, double rasterUnitsPerPixelX, double rasterUnitsPerPixelY, const QString &attributePrefix, int rasterBand, Qgis::ZonalStatistics stats )
45 : mRasterInterface( rasterInterface )
46 , mRasterCrs( rasterCrs )
47 , mCellSizeX( std::fabs( rasterUnitsPerPixelX ) )
48 , mCellSizeY( std::fabs( rasterUnitsPerPixelY ) )
49 , mRasterBand( rasterBand )
50 , mPolygonLayer( polygonLayer )
51 , mAttributePrefix( attributePrefix )
52 , mStatistics( stats )
53{
54}
55
57{
58 if ( !mRasterInterface )
59 {
61 }
62
63 if ( mRasterInterface->bandCount() < mRasterBand )
64 {
66 }
67
68 if ( !mPolygonLayer || mPolygonLayer->geometryType() != Qgis::GeometryType::Polygon )
69 {
71 }
72
73 QgsVectorDataProvider *vectorProvider = mPolygonLayer->dataProvider();
74 if ( !vectorProvider )
75 {
77 }
78
79 QMap<Qgis::ZonalStatistic, int> statFieldIndexes;
80
81 //add the new fields to the provider
82 int oldFieldCount = vectorProvider->fields().count();
83 QList<QgsField> newFieldList;
84 for ( Qgis::ZonalStatistic stat :
85 {
98 } )
99 {
100 if ( mStatistics & stat )
101 {
102 QString fieldName = getUniqueFieldName( mAttributePrefix + QgsZonalStatistics::shortName( stat ), newFieldList );
103 QgsField field( fieldName, QMetaType::Type::Double, QStringLiteral( "double precision" ) );
104 newFieldList.push_back( field );
105 statFieldIndexes.insert( stat, oldFieldCount + newFieldList.count() - 1 );
106 }
107 }
108
109 vectorProvider->addAttributes( newFieldList );
110
111 long featureCount = vectorProvider->featureCount();
112
113 QgsFeatureRequest request;
114 request.setNoAttributes();
115
116 request.setDestinationCrs( mRasterCrs, QgsProject::instance()->transformContext() );
117 QgsFeatureIterator fi = vectorProvider->getFeatures( request );
118 QgsFeature feature;
119
120 int featureCounter = 0;
121
122 QgsChangedAttributesMap changeMap;
123 while ( fi.nextFeature( feature ) )
124 {
125 ++featureCounter;
126 if ( feedback && feedback->isCanceled() )
127 {
128 break;
129 }
130
131 if ( feedback )
132 {
133 feedback->setProgress( 100.0 * static_cast< double >( featureCounter ) / featureCount );
134 }
135
136 QgsGeometry featureGeometry = feature.geometry();
137
138 QMap<Qgis::ZonalStatistic, QVariant> results = calculateStatistics( mRasterInterface, featureGeometry, mCellSizeX, mCellSizeY, mRasterBand, mStatistics );
139
140 if ( results.empty() )
141 continue;
142
143 QgsAttributeMap changeAttributeMap;
144 for ( const auto &result : results.toStdMap() )
145 {
146 changeAttributeMap.insert( statFieldIndexes.value( result.first ), result.second );
147 }
148
149 changeMap.insert( feature.id(), changeAttributeMap );
150 }
151
152 vectorProvider->changeAttributeValues( changeMap );
153 mPolygonLayer->updateFields();
154
155 if ( feedback )
156 {
157 if ( feedback->isCanceled() )
159
160 feedback->setProgress( 100 );
161 }
162
164}
165
166QString QgsZonalStatistics::getUniqueFieldName( const QString &fieldName, const QList<QgsField> &newFields )
167{
168 QgsVectorDataProvider *dp = mPolygonLayer->dataProvider();
169
170 if ( !dp->storageType().contains( QLatin1String( "ESRI Shapefile" ) ) )
171 {
172 return fieldName;
173 }
174
175 QList<QgsField> allFields = dp->fields().toList();
176 allFields.append( newFields );
177 QString shortName = fieldName.mid( 0, 10 );
178
179 bool found = false;
180 for ( const QgsField &field : std::as_const( allFields ) )
181 {
182 if ( shortName == field.name() )
183 {
184 found = true;
185 break;
186 }
187 }
188
189 if ( !found )
190 {
191 return shortName;
192 }
193
194 int n = 1;
195 shortName = QStringLiteral( "%1_%2" ).arg( fieldName.mid( 0, 8 ) ).arg( n );
196 found = true;
197 while ( found )
198 {
199 found = false;
200 for ( const QgsField &field : std::as_const( allFields ) )
201 {
202 if ( shortName == field.name() )
203 {
204 n += 1;
205 if ( n < 9 )
206 {
207 shortName = QStringLiteral( "%1_%2" ).arg( fieldName.mid( 0, 8 ) ).arg( n );
208 }
209 else
210 {
211 shortName = QStringLiteral( "%1_%2" ).arg( fieldName.mid( 0, 7 ) ).arg( n );
212 }
213 found = true;
214 }
215 }
216 }
217 return shortName;
218}
219
221{
222 switch ( statistic )
223 {
225 return QObject::tr( "Count" );
227 return QObject::tr( "Sum" );
229 return QObject::tr( "Mean" );
231 return QObject::tr( "Median" );
233 return QObject::tr( "St dev" );
235 return QObject::tr( "Minimum" );
237 return QObject::tr( "Maximum" );
239 return QObject::tr( "Minimum point" );
241 return QObject::tr( "Maximum point" );
243 return QObject::tr( "Range" );
245 return QObject::tr( "Minority" );
247 return QObject::tr( "Majority" );
249 return QObject::tr( "Variety" );
251 return QObject::tr( "Variance" );
254 return QString();
255 }
256 return QString();
257}
258
260{
261 switch ( statistic )
262 {
264 return QStringLiteral( "count" );
266 return QStringLiteral( "sum" );
268 return QStringLiteral( "mean" );
270 return QStringLiteral( "median" );
272 return QStringLiteral( "stdev" );
274 return QStringLiteral( "min" );
276 return QStringLiteral( "max" );
278 return QStringLiteral( "minpoint" );
280 return QStringLiteral( "maxpoint" );
282 return QStringLiteral( "range" );
284 return QStringLiteral( "minority" );
286 return QStringLiteral( "majority" );
288 return QStringLiteral( "variety" );
290 return QStringLiteral( "variance" );
293 return QString();
294 }
295 return QString();
296}
297
299QMap<int, QVariant> QgsZonalStatistics::calculateStatisticsInt( QgsRasterInterface *rasterInterface, const QgsGeometry &geometry, double cellSizeX, double cellSizeY, int rasterBand, Qgis::ZonalStatistics statistics )
300{
301 const auto result { QgsZonalStatistics::calculateStatistics( rasterInterface, geometry, cellSizeX, cellSizeY, rasterBand, statistics ) };
302 QMap<int, QVariant> pyResult;
303 for ( auto it = result.constBegin(); it != result.constEnd(); ++it )
304 {
305 pyResult.insert( static_cast< int >( it.key() ), it.value() );
306 }
307 return pyResult;
308}
310
311QMap<Qgis::ZonalStatistic, QVariant> QgsZonalStatistics::calculateStatistics( QgsRasterInterface *rasterInterface, const QgsGeometry &geometry, double cellSizeX, double cellSizeY, int rasterBand, Qgis::ZonalStatistics statistics )
312{
313 QMap<Qgis::ZonalStatistic, QVariant> results;
314
315 if ( geometry.isEmpty() )
316 return results;
317
318 QgsRectangle rasterBBox = rasterInterface->extent();
319
320 QgsRectangle featureRect = geometry.boundingBox().intersect( rasterBBox );
321
322 if ( featureRect.isEmpty() )
323 return results;
324
325 bool statsStoreValues = ( statistics & Qgis::ZonalStatistic::Median ) ||
326 ( statistics & Qgis::ZonalStatistic::StDev ) ||
327 ( statistics & Qgis::ZonalStatistic::Variance );
328 bool statsStoreValueCount = ( statistics & Qgis::ZonalStatistic::Minority ) ||
329 ( statistics & Qgis::ZonalStatistic::Majority );
330
331 FeatureStats featureStats( statsStoreValues, statsStoreValueCount );
332
333 int nCellsXProvider = rasterInterface->xSize();
334 int nCellsYProvider = rasterInterface->ySize();
335
336 int nCellsX, nCellsY;
337 QgsRectangle rasterBlockExtent;
338 QgsRasterAnalysisUtils::cellInfoForBBox( rasterBBox, featureRect, cellSizeX, cellSizeY, nCellsX, nCellsY, nCellsXProvider, nCellsYProvider, rasterBlockExtent );
339
340 featureStats.reset();
341 QgsRasterAnalysisUtils::statisticsFromMiddlePointTest( rasterInterface, rasterBand, geometry, nCellsX, nCellsY, cellSizeX, cellSizeY, rasterBlockExtent, [ &featureStats ]( double value, const QgsPointXY & point ) { featureStats.addValue( value, point ); } );
342
343 if ( featureStats.count <= 1 )
344 {
345 //the cell resolution is probably larger than the polygon area. We switch to precise pixel - polygon intersection in this case
346 featureStats.reset();
347 QgsRasterAnalysisUtils::statisticsFromPreciseIntersection( rasterInterface, rasterBand, geometry, nCellsX, nCellsY, cellSizeX, cellSizeY, rasterBlockExtent, [ &featureStats ]( double value, double weight, const QgsPointXY & point ) { featureStats.addValue( value, point, weight ); } );
348 }
349
350 // calculate the statistics
351
352 if ( statistics & Qgis::ZonalStatistic::Count )
353 results.insert( Qgis::ZonalStatistic::Count, QVariant( featureStats.count ) );
354 if ( statistics & Qgis::ZonalStatistic::Sum )
355 results.insert( Qgis::ZonalStatistic::Sum, QVariant( featureStats.sum ) );
356 if ( featureStats.count > 0 )
357 {
358 double mean = featureStats.sum / featureStats.count;
359 if ( statistics & Qgis::ZonalStatistic::Mean )
360 results.insert( Qgis::ZonalStatistic::Mean, QVariant( mean ) );
361 if ( statistics & Qgis::ZonalStatistic::Median )
362 {
363 std::sort( featureStats.values.begin(), featureStats.values.end() );
364 int size = featureStats.values.count();
365 bool even = ( size % 2 ) < 1;
366 double medianValue;
367 if ( even )
368 {
369 medianValue = ( featureStats.values.at( size / 2 - 1 ) + featureStats.values.at( size / 2 ) ) / 2;
370 }
371 else //odd
372 {
373 medianValue = featureStats.values.at( ( size + 1 ) / 2 - 1 );
374 }
375 results.insert( Qgis::ZonalStatistic::Median, QVariant( medianValue ) );
376 }
377 if ( statistics & Qgis::ZonalStatistic::StDev || statistics & Qgis::ZonalStatistic::Variance )
378 {
379 double sumSquared = 0;
380 for ( int i = 0; i < featureStats.values.count(); ++i )
381 {
382 double diff = featureStats.values.at( i ) - mean;
383 sumSquared += diff * diff;
384 }
385 double variance = sumSquared / featureStats.values.count();
386 if ( statistics & Qgis::ZonalStatistic::StDev )
387 {
388 double stdev = std::pow( variance, 0.5 );
389 results.insert( Qgis::ZonalStatistic::StDev, QVariant( stdev ) );
390 }
391 if ( statistics & Qgis::ZonalStatistic::Variance )
392 results.insert( Qgis::ZonalStatistic::Variance, QVariant( variance ) );
393 }
394 if ( statistics & Qgis::ZonalStatistic::Min )
395 results.insert( Qgis::ZonalStatistic::Min, QVariant( featureStats.min ) );
396 if ( statistics & Qgis::ZonalStatistic::Max )
397 results.insert( Qgis::ZonalStatistic::Max, QVariant( featureStats.max ) );
398 if ( statistics & Qgis::ZonalStatistic::MinimumPoint )
399 results.insert( Qgis::ZonalStatistic::MinimumPoint, QVariant( featureStats.minPoint ) );
400 if ( statistics & Qgis::ZonalStatistic::MaximumPoint )
401 results.insert( Qgis::ZonalStatistic::MaximumPoint, QVariant( featureStats.maxPoint ) );
402 if ( statistics & Qgis::ZonalStatistic::Range )
403 results.insert( Qgis::ZonalStatistic::Range, QVariant( featureStats.max - featureStats.min ) );
404 if ( statistics & Qgis::ZonalStatistic::Minority || statistics & Qgis::ZonalStatistic::Majority )
405 {
406 QList<int> vals = featureStats.valueCount.values();
407 std::sort( vals.begin(), vals.end() );
408 if ( statistics & Qgis::ZonalStatistic::Minority )
409 {
410 double minorityKey = featureStats.valueCount.key( vals.first() );
411 results.insert( Qgis::ZonalStatistic::Minority, QVariant( minorityKey ) );
412 }
413 if ( statistics & Qgis::ZonalStatistic::Majority )
414 {
415 double majKey = featureStats.valueCount.key( vals.last() );
416 results.insert( Qgis::ZonalStatistic::Majority, QVariant( majKey ) );
417 }
418 }
419 if ( statistics & Qgis::ZonalStatistic::Variety )
420 results.insert( Qgis::ZonalStatistic::Variety, QVariant( featureStats.valueCount.count() ) );
421 }
422
423 return results;
424}
ZonalStatistic
Statistics to be calculated during a zonal statistics operation.
Definition qgis.h:5406
@ StDev
Standard deviation of pixel values.
@ Mean
Mean of pixel values.
@ Median
Median of pixel values.
@ Max
Max of pixel values.
@ Variance
Variance of pixel values.
@ MinimumPoint
Pixel centroid for minimum pixel value.
@ Min
Min of pixel values.
@ Default
Default statistics.
@ Range
Range of pixel values (max - min)
@ Sum
Sum of pixel values.
@ Minority
Minority of pixel values.
@ All
All statistics. For QGIS 3.x this includes ONLY numeric statistics, but for 4.0 this will be extended...
@ Majority
Majority of pixel values.
@ Variety
Variety (count of distinct) pixel values.
@ MaximumPoint
Pixel centroid for maximum pixel value.
@ Count
Pixel count.
@ Polygon
Polygons.
ZonalStatisticResult
Zonal statistics result codes.
Definition qgis.h:5441
@ Canceled
Algorithm was canceled.
@ LayerTypeWrong
Layer is not a polygon layer.
@ RasterBandInvalid
The raster band does not exist on the raster layer.
@ RasterInvalid
Raster layer is invalid.
@ LayerInvalid
Layer is invalid.
QFlags< ZonalStatistic > ZonalStatistics
Statistics to be calculated during a zonal statistics operation.
Definition qgis.h:5432
This class represents a coordinate reference system (CRS).
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.
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition qgsfeature.h:58
QgsFeatureId id
Definition qgsfeature.h:66
QgsGeometry geometry
Definition qgsfeature.h:69
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition qgsfeedback.h:44
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
int count
Definition qgsfields.h:50
QList< QgsField > toList() const
Utility function to return a list of QgsField instances.
A geometry is the spatial representation of a feature.
bool isEmpty() const
Returns true if the geometry is empty (eg a linestring with no vertices, or a collection with no geom...
QgsRectangle boundingBox() const
Returns the bounding box of the geometry.
A class to represent a 2D point.
Definition qgspointxy.h:60
static QgsProject * instance()
Returns the QgsProject singleton instance.
Base class for processing filters like renderers, reprojector, resampler etc.
virtual int xSize() const
Gets raster size.
virtual int bandCount() const =0
Gets number of bands.
virtual int ySize() const
virtual QgsRectangle extent() const
Gets the extent of the interface.
Represents a raster layer.
A rectangle specified with double values.
QgsRectangle intersect(const QgsRectangle &rect) const
Returns the intersection with the given rectangle.
This is the base class for vector data providers.
long long featureCount() const override=0
Number of features in the layer.
virtual QString storageType() const
Returns the permanent storage type for this layer as a friendly name.
virtual bool changeAttributeValues(const QgsChangedAttributesMap &attr_map)
Changes attribute values of existing features.
QgsFields fields() const override=0
Returns the fields associated with this data provider.
QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest()) const override=0
Query the provider for features specified in request.
virtual bool addAttributes(const QList< QgsField > &attributes)
Adds new attributes to the provider.
Represents a vector layer which manages a vector based data sets.
void updateFields()
Will regenerate the fields property of this layer by obtaining all fields from the dataProvider,...
QgsVectorDataProvider * dataProvider() FINAL
Returns the layer's data provider, it may be nullptr.
Q_INVOKABLE Qgis::GeometryType geometryType() const
Returns point, line or polygon.
A class that calculates raster statistics (count, sum, mean) for a polygon or multipolygon layer and ...
QgsZonalStatistics(QgsVectorLayer *polygonLayer, QgsRasterLayer *rasterLayer, const QString &attributePrefix=QString(), int rasterBand=1, Qgis::ZonalStatistics stats=Qgis::ZonalStatistic::Default)
Convenience constructor for QgsZonalStatistics, using an input raster layer.
Qgis::ZonalStatisticResult calculateStatistics(QgsFeedback *feedback)
Runs the calculation.
static QString displayName(Qgis::ZonalStatistic statistic)
Returns the friendly display name for a statistic.
static QString shortName(Qgis::ZonalStatistic statistic)
Returns a short, friendly display name for a statistic, suitable for use in a field name.
QMap< int, QVariant > QgsAttributeMap
QMap< QgsFeatureId, QgsAttributeMap > QgsChangedAttributesMap
const QgsCoordinateReferenceSystem & crs