QGIS API Documentation  3.21.0-Master (5b68dc587e)
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"
23 #include "qgsvectordataprovider.h"
24 #include "qgsvectorlayer.h"
26 #include "qgsrasterdataprovider.h"
27 #include "qgsrasterlayer.h"
28 #include "qgslogger.h"
29 #include "qgsproject.h"
30 
31 #include <QFile>
32 
33 QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer *polygonLayer, QgsRasterLayer *rasterLayer, const QString &attributePrefix, int rasterBand, QgsZonalStatistics::Statistics stats )
34  : QgsZonalStatistics( polygonLayer,
35  rasterLayer ? rasterLayer->dataProvider() : nullptr,
36  rasterLayer ? rasterLayer->crs() : QgsCoordinateReferenceSystem(),
37  rasterLayer ? rasterLayer->rasterUnitsPerPixelX() : 0,
38  rasterLayer ? rasterLayer->rasterUnitsPerPixelY() : 0,
39  attributePrefix,
40  rasterBand,
41  stats )
42 {
43 }
44 
46  const QgsCoordinateReferenceSystem &rasterCrs, double rasterUnitsPerPixelX, double rasterUnitsPerPixelY, const QString &attributePrefix, int rasterBand, QgsZonalStatistics::Statistics stats )
47  : mRasterInterface( rasterInterface )
48  , mRasterCrs( rasterCrs )
49  , mCellSizeX( std::fabs( rasterUnitsPerPixelX ) )
50  , mCellSizeY( std::fabs( rasterUnitsPerPixelY ) )
51  , mRasterBand( rasterBand )
52  , mPolygonLayer( polygonLayer )
53  , mAttributePrefix( attributePrefix )
54  , mStatistics( stats )
55 {
56 }
57 
59 {
60  if ( !mRasterInterface )
61  {
62  return RasterInvalid;
63  }
64 
65  if ( mRasterInterface->bandCount() < mRasterBand )
66  {
67  return RasterBandInvalid;
68  }
69 
70  if ( !mPolygonLayer || mPolygonLayer->geometryType() != QgsWkbTypes::PolygonGeometry )
71  {
72  return LayerTypeWrong;
73  }
74 
75  QgsVectorDataProvider *vectorProvider = mPolygonLayer->dataProvider();
76  if ( !vectorProvider )
77  {
78  return LayerInvalid;
79  }
80 
81  QMap<QgsZonalStatistics::Statistic, int> statFieldIndexes;
82 
83  //add the new fields to the provider
84  int oldFieldCount = vectorProvider->fields().count();
85  QList<QgsField> newFieldList;
87  {
100  } )
101  {
102  if ( mStatistics & stat )
103  {
104  QString fieldName = getUniqueFieldName( mAttributePrefix + QgsZonalStatistics::shortName( stat ), newFieldList );
105  QgsField field( fieldName, QVariant::Double, QStringLiteral( "double precision" ) );
106  newFieldList.push_back( field );
107  statFieldIndexes.insert( stat, oldFieldCount + newFieldList.count() - 1 );
108  }
109  }
110 
111  vectorProvider->addAttributes( newFieldList );
112 
113  long featureCount = vectorProvider->featureCount();
114 
115  QgsFeatureRequest request;
116  request.setNoAttributes();
117 
118  request.setDestinationCrs( mRasterCrs, QgsProject::instance()->transformContext() );
119  QgsFeatureIterator fi = vectorProvider->getFeatures( request );
120  QgsFeature feature;
121 
122  int featureCounter = 0;
123 
124  QgsChangedAttributesMap changeMap;
125  while ( fi.nextFeature( feature ) )
126  {
127  ++featureCounter;
128  if ( feedback && feedback->isCanceled() )
129  {
130  break;
131  }
132 
133  if ( feedback )
134  {
135  feedback->setProgress( 100.0 * static_cast< double >( featureCounter ) / featureCount );
136  }
137 
138  QgsGeometry featureGeometry = feature.geometry();
139 
140  QMap<QgsZonalStatistics::Statistic, QVariant> results = calculateStatistics( mRasterInterface, featureGeometry, mCellSizeX, mCellSizeY, mRasterBand, mStatistics );
141 
142  if ( results.empty() )
143  continue;
144 
145  QgsAttributeMap changeAttributeMap;
146  for ( const auto &result : results.toStdMap() )
147  {
148  changeAttributeMap.insert( statFieldIndexes.value( result.first ), result.second );
149  }
150 
151  changeMap.insert( feature.id(), changeAttributeMap );
152  }
153 
154  vectorProvider->changeAttributeValues( changeMap );
155  mPolygonLayer->updateFields();
156 
157  if ( feedback )
158  {
159  if ( feedback->isCanceled() )
160  return Canceled;
161 
162  feedback->setProgress( 100 );
163  }
164 
165  return Success;
166 }
167 
168 QString QgsZonalStatistics::getUniqueFieldName( const QString &fieldName, const QList<QgsField> &newFields )
169 {
170  QgsVectorDataProvider *dp = mPolygonLayer->dataProvider();
171 
172  if ( !dp->storageType().contains( QLatin1String( "ESRI Shapefile" ) ) )
173  {
174  return fieldName;
175  }
176 
177  QList<QgsField> allFields = dp->fields().toList();
178  allFields.append( newFields );
179  QString shortName = fieldName.mid( 0, 10 );
180 
181  bool found = false;
182  for ( const QgsField &field : std::as_const( allFields ) )
183  {
184  if ( shortName == field.name() )
185  {
186  found = true;
187  break;
188  }
189  }
190 
191  if ( !found )
192  {
193  return shortName;
194  }
195 
196  int n = 1;
197  shortName = QStringLiteral( "%1_%2" ).arg( fieldName.mid( 0, 8 ) ).arg( n );
198  found = true;
199  while ( found )
200  {
201  found = false;
202  for ( const QgsField &field : std::as_const( allFields ) )
203  {
204  if ( shortName == field.name() )
205  {
206  n += 1;
207  if ( n < 9 )
208  {
209  shortName = QStringLiteral( "%1_%2" ).arg( fieldName.mid( 0, 8 ) ).arg( n );
210  }
211  else
212  {
213  shortName = QStringLiteral( "%1_%2" ).arg( fieldName.mid( 0, 7 ) ).arg( n );
214  }
215  found = true;
216  }
217  }
218  }
219  return shortName;
220 }
221 
223 {
224  switch ( statistic )
225  {
226  case Count:
227  return QObject::tr( "Count" );
228  case Sum:
229  return QObject::tr( "Sum" );
230  case Mean:
231  return QObject::tr( "Mean" );
232  case Median:
233  return QObject::tr( "Median" );
234  case StDev:
235  return QObject::tr( "St dev" );
236  case Min:
237  return QObject::tr( "Minimum" );
238  case Max:
239  return QObject::tr( "Maximum" );
240  case Range:
241  return QObject::tr( "Range" );
242  case Minority:
243  return QObject::tr( "Minority" );
244  case Majority:
245  return QObject::tr( "Majority" );
246  case Variety:
247  return QObject::tr( "Variety" );
248  case Variance:
249  return QObject::tr( "Variance" );
250  case All:
251  return QString();
252  }
253  return QString();
254 }
255 
257 {
258  switch ( statistic )
259  {
260  case Count:
261  return QStringLiteral( "count" );
262  case Sum:
263  return QStringLiteral( "sum" );
264  case Mean:
265  return QStringLiteral( "mean" );
266  case Median:
267  return QStringLiteral( "median" );
268  case StDev:
269  return QStringLiteral( "stdev" );
270  case Min:
271  return QStringLiteral( "min" );
272  case Max:
273  return QStringLiteral( "max" );
274  case Range:
275  return QStringLiteral( "range" );
276  case Minority:
277  return QStringLiteral( "minority" );
278  case Majority:
279  return QStringLiteral( "majority" );
280  case Variety:
281  return QStringLiteral( "variety" );
282  case Variance:
283  return QStringLiteral( "variance" );
284  case All:
285  return QString();
286  }
287  return QString();
288 }
289 
291 QMap<int, QVariant> QgsZonalStatistics::calculateStatisticsInt( QgsRasterInterface *rasterInterface, const QgsGeometry &geometry, double cellSizeX, double cellSizeY, int rasterBand, QgsZonalStatistics::Statistics statistics )
292 {
293  const auto result { QgsZonalStatistics::calculateStatistics( rasterInterface, geometry, cellSizeX, cellSizeY, rasterBand, statistics ) };
294  QMap<int, QVariant> pyResult;
295  for ( auto it = result.constBegin(); it != result.constEnd(); ++it )
296  {
297  pyResult.insert( it.key(), it.value() );
298  }
299  return pyResult;
300 }
302 
303 QMap<QgsZonalStatistics::Statistic, QVariant> QgsZonalStatistics::calculateStatistics( QgsRasterInterface *rasterInterface, const QgsGeometry &geometry, double cellSizeX, double cellSizeY, int rasterBand, QgsZonalStatistics::Statistics statistics )
304 {
305  QMap<QgsZonalStatistics::Statistic, QVariant> results;
306 
307  if ( geometry.isEmpty() )
308  return results;
309 
310  QgsRectangle rasterBBox = rasterInterface->extent();
311 
312  QgsRectangle featureRect = geometry.boundingBox().intersect( rasterBBox );
313 
314  if ( featureRect.isEmpty() )
315  return results;
316 
317  bool statsStoreValues = ( statistics & QgsZonalStatistics::Median ) ||
318  ( statistics & QgsZonalStatistics::StDev ) ||
319  ( statistics & QgsZonalStatistics::Variance );
320  bool statsStoreValueCount = ( statistics & QgsZonalStatistics::Minority ) ||
321  ( statistics & QgsZonalStatistics::Majority );
322 
323  FeatureStats featureStats( statsStoreValues, statsStoreValueCount );
324 
325  int nCellsXProvider = rasterInterface->xSize();
326  int nCellsYProvider = rasterInterface->ySize();
327 
328  int nCellsX, nCellsY;
329  QgsRectangle rasterBlockExtent;
330  QgsRasterAnalysisUtils::cellInfoForBBox( rasterBBox, featureRect, cellSizeX, cellSizeY, nCellsX, nCellsY, nCellsXProvider, nCellsYProvider, rasterBlockExtent );
331 
332  featureStats.reset();
333  QgsRasterAnalysisUtils::statisticsFromMiddlePointTest( rasterInterface, rasterBand, geometry, nCellsX, nCellsY, cellSizeX, cellSizeY, rasterBlockExtent, [ &featureStats ]( double value ) { featureStats.addValue( value ); } );
334 
335  if ( featureStats.count <= 1 )
336  {
337  //the cell resolution is probably larger than the polygon area. We switch to precise pixel - polygon intersection in this case
338  featureStats.reset();
339  QgsRasterAnalysisUtils::statisticsFromPreciseIntersection( rasterInterface, rasterBand, geometry, nCellsX, nCellsY, cellSizeX, cellSizeY, rasterBlockExtent, [ &featureStats ]( double value, double weight ) { featureStats.addValue( value, weight ); } );
340  }
341 
342  // calculate the statistics
343 
344  if ( statistics & QgsZonalStatistics::Count )
345  results.insert( QgsZonalStatistics::Count, QVariant( featureStats.count ) );
346  if ( statistics & QgsZonalStatistics::Sum )
347  results.insert( QgsZonalStatistics::Sum, QVariant( featureStats.sum ) );
348  if ( featureStats.count > 0 )
349  {
350  double mean = featureStats.sum / featureStats.count;
351  if ( statistics & QgsZonalStatistics::Mean )
352  results.insert( QgsZonalStatistics::Mean, QVariant( mean ) );
353  if ( statistics & QgsZonalStatistics::Median )
354  {
355  std::sort( featureStats.values.begin(), featureStats.values.end() );
356  int size = featureStats.values.count();
357  bool even = ( size % 2 ) < 1;
358  double medianValue;
359  if ( even )
360  {
361  medianValue = ( featureStats.values.at( size / 2 - 1 ) + featureStats.values.at( size / 2 ) ) / 2;
362  }
363  else //odd
364  {
365  medianValue = featureStats.values.at( ( size + 1 ) / 2 - 1 );
366  }
367  results.insert( QgsZonalStatistics::Median, QVariant( medianValue ) );
368  }
369  if ( statistics & QgsZonalStatistics::StDev || statistics & QgsZonalStatistics::Variance )
370  {
371  double sumSquared = 0;
372  for ( int i = 0; i < featureStats.values.count(); ++i )
373  {
374  double diff = featureStats.values.at( i ) - mean;
375  sumSquared += diff * diff;
376  }
377  double variance = sumSquared / featureStats.values.count();
378  if ( statistics & QgsZonalStatistics::StDev )
379  {
380  double stdev = std::pow( variance, 0.5 );
381  results.insert( QgsZonalStatistics::StDev, QVariant( stdev ) );
382  }
383  if ( statistics & QgsZonalStatistics::Variance )
384  results.insert( QgsZonalStatistics::Variance, QVariant( variance ) );
385  }
386  if ( statistics & QgsZonalStatistics::Min )
387  results.insert( QgsZonalStatistics::Min, QVariant( featureStats.min ) );
388  if ( statistics & QgsZonalStatistics::Max )
389  results.insert( QgsZonalStatistics::Max, QVariant( featureStats.max ) );
390  if ( statistics & QgsZonalStatistics::Range )
391  results.insert( QgsZonalStatistics::Range, QVariant( featureStats.max - featureStats.min ) );
392  if ( statistics & QgsZonalStatistics::Minority || statistics & QgsZonalStatistics::Majority )
393  {
394  QList<int> vals = featureStats.valueCount.values();
395  std::sort( vals.begin(), vals.end() );
396  if ( statistics & QgsZonalStatistics::Minority )
397  {
398  double minorityKey = featureStats.valueCount.key( vals.first() );
399  results.insert( QgsZonalStatistics::Minority, QVariant( minorityKey ) );
400  }
401  if ( statistics & QgsZonalStatistics::Majority )
402  {
403  double majKey = featureStats.valueCount.key( vals.last() );
404  results.insert( QgsZonalStatistics::Majority, QVariant( majKey ) );
405  }
406  }
407  if ( statistics & QgsZonalStatistics::Variety )
408  results.insert( QgsZonalStatistics::Variety, QVariant( featureStats.valueCount.count() ) );
409  }
410 
411  return results;
412 }
This class represents a coordinate reference system (CRS).
Wrapper for iterator of features from vector data provider or vector layer.
bool nextFeature(QgsFeature &f)
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:56
QgsGeometry geometry
Definition: qgsfeature.h:67
Q_GADGET QgsFeatureId id
Definition: qgsfeature.h:64
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition: qgsfeedback.h:45
bool isCanceled() const SIP_HOLDGIL
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
Encapsulate a field in an attribute table or data source.
Definition: qgsfield.h:51
QString name
Definition: qgsfield.h:60
QList< QgsField > toList() const
Utility function to return a list of QgsField instances.
Definition: qgsfields.cpp:212
int count() const
Returns number of items.
Definition: qgsfields.cpp:133
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:124
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.
static QgsProject * instance()
Returns the QgsProject singleton instance.
Definition: qgsproject.cpp:467
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.
Definition: qgsrectangle.h:42
bool isEmpty() const
Returns true if the rectangle is empty.
Definition: qgsrectangle.h:469
QgsRectangle intersect(const QgsRectangle &rect) const
Returns the intersection with the given rectangle.
Definition: qgsrectangle.h:333
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.
Q_INVOKABLE QgsWkbTypes::GeometryType geometryType() const
Returns point, line or polygon.
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.
A class that calculates raster statistics (count, sum, mean) for a polygon or multipolygon layer and ...
static QString shortName(QgsZonalStatistics::Statistic statistic)
Returns a short, friendly display name for a statistic, suitable for use in a field name.
Statistic
Enumeration of flags that specify statistics to be calculated.
@ Minority
Minority of pixel values.
@ Variety
Variety (count of distinct) pixel values.
@ Variance
Variance of pixel values.
@ Mean
Mean of pixel values.
@ Majority
Majority of pixel values.
@ Range
Range of pixel values (max - min)
@ Min
Min of pixel values.
@ Sum
Sum of pixel values.
@ StDev
Standard deviation of pixel values.
@ Max
Max of pixel values.
@ Median
Median of pixel values.
Result
Error codes for calculation.
@ RasterInvalid
Raster layer is invalid.
@ RasterBandInvalid
The raster band does not exist on the raster layer.
@ LayerTypeWrong
Layer is not a polygon layer.
@ Canceled
Algorithm was canceled.
@ LayerInvalid
Layer is invalid.
QgsZonalStatistics::Result calculateStatistics(QgsFeedback *feedback)
Runs the calculation.
static QString displayName(QgsZonalStatistics::Statistic statistic)
Returns the friendly display name for a statistic.
QgsZonalStatistics(QgsVectorLayer *polygonLayer, QgsRasterLayer *rasterLayer, const QString &attributePrefix=QString(), int rasterBand=1, QgsZonalStatistics::Statistics stats=QgsZonalStatistics::Statistics(QgsZonalStatistics::Count|QgsZonalStatistics::Sum|QgsZonalStatistics::Mean))
Convenience constructor for QgsZonalStatistics, using an input raster layer.
QMap< int, QVariant > QgsAttributeMap
Definition: qgsattributes.h:38
QMap< QgsFeatureId, QgsAttributeMap > QgsChangedAttributesMap
Definition: qgsfeature.h:820
const QgsField & field
Definition: qgsfield.h:463
const QgsCoordinateReferenceSystem & crs