QGIS API Documentation  3.9.0-Master (d9ef585e47)
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 ( !mPolygonLayer || mPolygonLayer->geometryType() != QgsWkbTypes::PolygonGeometry )
61  {
62  return 1;
63  }
64 
65  QgsVectorDataProvider *vectorProvider = mPolygonLayer->dataProvider();
66  if ( !vectorProvider )
67  {
68  return 2;
69  }
70 
71  if ( !mRasterInterface )
72  {
73  return 3;
74  }
75 
76  if ( mRasterInterface->bandCount() < mRasterBand )
77  {
78  return 4;
79  }
80 
81  //get geometry info about raster layer
82  int nCellsXProvider = mRasterInterface->xSize();
83  int nCellsYProvider = mRasterInterface->ySize();
84 
85  QgsRectangle rasterBBox = mRasterInterface->extent();
86 
87  //add the new fields to the provider
88  QList<QgsField> newFieldList;
89  QString countFieldName;
90  if ( mStatistics & QgsZonalStatistics::Count )
91  {
92  countFieldName = getUniqueFieldName( mAttributePrefix + "count", newFieldList );
93  QgsField countField( countFieldName, QVariant::Double, QStringLiteral( "double precision" ) );
94  newFieldList.push_back( countField );
95  }
96  QString sumFieldName;
97  if ( mStatistics & QgsZonalStatistics::Sum )
98  {
99  sumFieldName = getUniqueFieldName( mAttributePrefix + "sum", newFieldList );
100  QgsField sumField( sumFieldName, QVariant::Double, QStringLiteral( "double precision" ) );
101  newFieldList.push_back( sumField );
102  }
103  QString meanFieldName;
104  if ( mStatistics & QgsZonalStatistics::Mean )
105  {
106  meanFieldName = getUniqueFieldName( mAttributePrefix + "mean", newFieldList );
107  QgsField meanField( meanFieldName, QVariant::Double, QStringLiteral( "double precision" ) );
108  newFieldList.push_back( meanField );
109  }
110  QString medianFieldName;
111  if ( mStatistics & QgsZonalStatistics::Median )
112  {
113  medianFieldName = getUniqueFieldName( mAttributePrefix + "median", newFieldList );
114  QgsField medianField( medianFieldName, QVariant::Double, QStringLiteral( "double precision" ) );
115  newFieldList.push_back( medianField );
116  }
117  QString stdevFieldName;
118  if ( mStatistics & QgsZonalStatistics::StDev )
119  {
120  stdevFieldName = getUniqueFieldName( mAttributePrefix + "stdev", newFieldList );
121  QgsField stdField( stdevFieldName, QVariant::Double, QStringLiteral( "double precision" ) );
122  newFieldList.push_back( stdField );
123  }
124  QString minFieldName;
125  if ( mStatistics & QgsZonalStatistics::Min )
126  {
127  minFieldName = getUniqueFieldName( mAttributePrefix + "min", newFieldList );
128  QgsField minField( minFieldName, QVariant::Double, QStringLiteral( "double precision" ) );
129  newFieldList.push_back( minField );
130  }
131  QString maxFieldName;
132  if ( mStatistics & QgsZonalStatistics::Max )
133  {
134  maxFieldName = getUniqueFieldName( mAttributePrefix + "max", newFieldList );
135  QgsField maxField( maxFieldName, QVariant::Double, QStringLiteral( "double precision" ) );
136  newFieldList.push_back( maxField );
137  }
138  QString rangeFieldName;
139  if ( mStatistics & QgsZonalStatistics::Range )
140  {
141  rangeFieldName = getUniqueFieldName( mAttributePrefix + "range", newFieldList );
142  QgsField rangeField( rangeFieldName, QVariant::Double, QStringLiteral( "double precision" ) );
143  newFieldList.push_back( rangeField );
144  }
145  QString minorityFieldName;
146  if ( mStatistics & QgsZonalStatistics::Minority )
147  {
148  minorityFieldName = getUniqueFieldName( mAttributePrefix + "minority", newFieldList );
149  QgsField minorityField( minorityFieldName, QVariant::Double, QStringLiteral( "double precision" ) );
150  newFieldList.push_back( minorityField );
151  }
152  QString majorityFieldName;
153  if ( mStatistics & QgsZonalStatistics::Majority )
154  {
155  majorityFieldName = getUniqueFieldName( mAttributePrefix + "majority", newFieldList );
156  QgsField majField( majorityFieldName, QVariant::Double, QStringLiteral( "double precision" ) );
157  newFieldList.push_back( majField );
158  }
159  QString varietyFieldName;
160  if ( mStatistics & QgsZonalStatistics::Variety )
161  {
162  varietyFieldName = getUniqueFieldName( mAttributePrefix + "variety", newFieldList );
163  QgsField varietyField( varietyFieldName, QVariant::Int, QStringLiteral( "int" ) );
164  newFieldList.push_back( varietyField );
165  }
166  QString varianceFieldName;
167  if ( mStatistics & QgsZonalStatistics::Variance )
168  {
169  varianceFieldName = getUniqueFieldName( mAttributePrefix + "variance", newFieldList );
170  QgsField varianceField( varianceFieldName, QVariant::Double, QStringLiteral( "double precision" ) );
171  newFieldList.push_back( varianceField );
172  }
173  vectorProvider->addAttributes( newFieldList );
174 
175  //index of the new fields
176  int countIndex = mStatistics & QgsZonalStatistics::Count ? vectorProvider->fieldNameIndex( countFieldName ) : -1;
177  int sumIndex = mStatistics & QgsZonalStatistics::Sum ? vectorProvider->fieldNameIndex( sumFieldName ) : -1;
178  int meanIndex = mStatistics & QgsZonalStatistics::Mean ? vectorProvider->fieldNameIndex( meanFieldName ) : -1;
179  int medianIndex = mStatistics & QgsZonalStatistics::Median ? vectorProvider->fieldNameIndex( medianFieldName ) : -1;
180  int stdevIndex = mStatistics & QgsZonalStatistics::StDev ? vectorProvider->fieldNameIndex( stdevFieldName ) : -1;
181  int minIndex = mStatistics & QgsZonalStatistics::Min ? vectorProvider->fieldNameIndex( minFieldName ) : -1;
182  int maxIndex = mStatistics & QgsZonalStatistics::Max ? vectorProvider->fieldNameIndex( maxFieldName ) : -1;
183  int rangeIndex = mStatistics & QgsZonalStatistics::Range ? vectorProvider->fieldNameIndex( rangeFieldName ) : -1;
184  int minorityIndex = mStatistics & QgsZonalStatistics::Minority ? vectorProvider->fieldNameIndex( minorityFieldName ) : -1;
185  int majorityIndex = mStatistics & QgsZonalStatistics::Majority ? vectorProvider->fieldNameIndex( majorityFieldName ) : -1;
186  int varietyIndex = mStatistics & QgsZonalStatistics::Variety ? vectorProvider->fieldNameIndex( varietyFieldName ) : -1;
187  int varianceIndex = mStatistics & QgsZonalStatistics::Variance ? vectorProvider->fieldNameIndex( varianceFieldName ) : -1;
188 
189  if ( ( mStatistics & QgsZonalStatistics::Count && countIndex == -1 )
190  || ( mStatistics & QgsZonalStatistics::Sum && sumIndex == -1 )
191  || ( mStatistics & QgsZonalStatistics::Mean && meanIndex == -1 )
192  || ( mStatistics & QgsZonalStatistics::Median && medianIndex == -1 )
193  || ( mStatistics & QgsZonalStatistics::StDev && stdevIndex == -1 )
194  || ( mStatistics & QgsZonalStatistics::Min && minIndex == -1 )
195  || ( mStatistics & QgsZonalStatistics::Max && maxIndex == -1 )
196  || ( mStatistics & QgsZonalStatistics::Range && rangeIndex == -1 )
197  || ( mStatistics & QgsZonalStatistics::Minority && minorityIndex == -1 )
198  || ( mStatistics & QgsZonalStatistics::Majority && majorityIndex == -1 )
199  || ( mStatistics & QgsZonalStatistics::Variety && varietyIndex == -1 )
200  || ( mStatistics & QgsZonalStatistics::Variance && varianceIndex == -1 )
201  )
202  {
203  //failed to create a required field
204  return 8;
205  }
206 
207  //progress dialog
208  long featureCount = vectorProvider->featureCount();
209 
210  //iterate over each polygon
211  QgsFeatureRequest request;
212  request.setNoAttributes();
213  request.setDestinationCrs( mRasterCrs, QgsProject::instance()->transformContext() );
214  QgsFeatureIterator fi = vectorProvider->getFeatures( request );
215  QgsFeature f;
216 
217  bool statsStoreValues = ( mStatistics & QgsZonalStatistics::Median ) ||
218  ( mStatistics & QgsZonalStatistics::StDev ) ||
219  ( mStatistics & QgsZonalStatistics::Variance );
220  bool statsStoreValueCount = ( mStatistics & QgsZonalStatistics::Minority ) ||
221  ( mStatistics & QgsZonalStatistics::Majority );
222 
223  FeatureStats featureStats( statsStoreValues, statsStoreValueCount );
224  int featureCounter = 0;
225 
226  QgsChangedAttributesMap changeMap;
227  while ( fi.nextFeature( f ) )
228  {
229  if ( feedback && feedback->isCanceled() )
230  {
231  break;
232  }
233 
234  if ( feedback )
235  {
236  feedback->setProgress( 100.0 * static_cast< double >( featureCounter ) / featureCount );
237  }
238 
239  if ( !f.hasGeometry() )
240  {
241  ++featureCounter;
242  continue;
243  }
244  QgsGeometry featureGeometry = f.geometry();
245 
246  QgsRectangle featureRect = featureGeometry.boundingBox().intersect( rasterBBox );
247  if ( featureRect.isEmpty() )
248  {
249  ++featureCounter;
250  continue;
251  }
252 
253  int nCellsX, nCellsY;
254  QgsRectangle rasterBlockExtent;
255  QgsRasterAnalysisUtils::cellInfoForBBox( rasterBBox, featureRect, mCellSizeX, mCellSizeY, nCellsX, nCellsY, nCellsXProvider, nCellsYProvider, rasterBlockExtent );
256 
257  featureStats.reset();
258  QgsRasterAnalysisUtils::statisticsFromMiddlePointTest( mRasterInterface, mRasterBand, featureGeometry, nCellsX, nCellsY, mCellSizeX, mCellSizeY,
259  rasterBlockExtent, [ &featureStats ]( double value ) { featureStats.addValue( value ); } );
260 
261  if ( featureStats.count <= 1 )
262  {
263  //the cell resolution is probably larger than the polygon area. We switch to precise pixel - polygon intersection in this case
264  featureStats.reset();
265  QgsRasterAnalysisUtils::statisticsFromPreciseIntersection( mRasterInterface, mRasterBand, featureGeometry, nCellsX, nCellsY, mCellSizeX, mCellSizeY,
266  rasterBlockExtent, [ &featureStats ]( double value, double weight ) { featureStats.addValue( value, weight ); } );
267  }
268 
269  //write the statistics value to the vector data provider
270  QgsAttributeMap changeAttributeMap;
271  if ( mStatistics & QgsZonalStatistics::Count )
272  changeAttributeMap.insert( countIndex, QVariant( featureStats.count ) );
273  if ( mStatistics & QgsZonalStatistics::Sum )
274  changeAttributeMap.insert( sumIndex, QVariant( featureStats.sum ) );
275  if ( featureStats.count > 0 )
276  {
277  double mean = featureStats.sum / featureStats.count;
278  if ( mStatistics & QgsZonalStatistics::Mean )
279  changeAttributeMap.insert( meanIndex, QVariant( mean ) );
280  if ( mStatistics & QgsZonalStatistics::Median )
281  {
282  std::sort( featureStats.values.begin(), featureStats.values.end() );
283  int size = featureStats.values.count();
284  bool even = ( size % 2 ) < 1;
285  double medianValue;
286  if ( even )
287  {
288  medianValue = ( featureStats.values.at( size / 2 - 1 ) + featureStats.values.at( size / 2 ) ) / 2;
289  }
290  else //odd
291  {
292  medianValue = featureStats.values.at( ( size + 1 ) / 2 - 1 );
293  }
294  changeAttributeMap.insert( medianIndex, QVariant( medianValue ) );
295  }
296  if ( mStatistics & QgsZonalStatistics::StDev || mStatistics & QgsZonalStatistics::Variance )
297  {
298  double sumSquared = 0;
299  for ( int i = 0; i < featureStats.values.count(); ++i )
300  {
301  double diff = featureStats.values.at( i ) - mean;
302  sumSquared += diff * diff;
303  }
304  double variance = sumSquared / featureStats.values.count();
305  if ( mStatistics & QgsZonalStatistics::StDev )
306  {
307  double stdev = std::pow( variance, 0.5 );
308  changeAttributeMap.insert( stdevIndex, QVariant( stdev ) );
309  }
310  if ( mStatistics & QgsZonalStatistics::Variance )
311  changeAttributeMap.insert( varianceIndex, QVariant( variance ) );
312  }
313  if ( mStatistics & QgsZonalStatistics::Min )
314  changeAttributeMap.insert( minIndex, QVariant( featureStats.min ) );
315  if ( mStatistics & QgsZonalStatistics::Max )
316  changeAttributeMap.insert( maxIndex, QVariant( featureStats.max ) );
317  if ( mStatistics & QgsZonalStatistics::Range )
318  changeAttributeMap.insert( rangeIndex, QVariant( featureStats.max - featureStats.min ) );
319  if ( mStatistics & QgsZonalStatistics::Minority || mStatistics & QgsZonalStatistics::Majority )
320  {
321  QList<int> vals = featureStats.valueCount.values();
322  std::sort( vals.begin(), vals.end() );
323  if ( mStatistics & QgsZonalStatistics::Minority )
324  {
325  double minorityKey = featureStats.valueCount.key( vals.first() );
326  changeAttributeMap.insert( minorityIndex, QVariant( minorityKey ) );
327  }
328  if ( mStatistics & QgsZonalStatistics::Majority )
329  {
330  double majKey = featureStats.valueCount.key( vals.last() );
331  changeAttributeMap.insert( majorityIndex, QVariant( majKey ) );
332  }
333  }
334  if ( mStatistics & QgsZonalStatistics::Variety )
335  changeAttributeMap.insert( varietyIndex, QVariant( featureStats.valueCount.count() ) );
336  }
337 
338  changeMap.insert( f.id(), changeAttributeMap );
339  ++featureCounter;
340  }
341 
342  vectorProvider->changeAttributeValues( changeMap );
343 
344  if ( feedback )
345  {
346  feedback->setProgress( 100 );
347  }
348 
349  mPolygonLayer->updateFields();
350 
351  if ( feedback && feedback->isCanceled() )
352  {
353  return 9;
354  }
355 
356  return 0;
357 }
358 
359 QString QgsZonalStatistics::getUniqueFieldName( const QString &fieldName, const QList<QgsField> &newFields )
360 {
361  QgsVectorDataProvider *dp = mPolygonLayer->dataProvider();
362 
363  if ( !dp->storageType().contains( QLatin1String( "ESRI Shapefile" ) ) )
364  {
365  return fieldName;
366  }
367 
368  QList<QgsField> allFields = dp->fields().toList();
369  allFields.append( newFields );
370  QString shortName = fieldName.mid( 0, 10 );
371 
372  bool found = false;
373  for ( int idx = 0; idx < allFields.count(); ++idx )
374  {
375  if ( shortName == allFields.at( idx ).name() )
376  {
377  found = true;
378  break;
379  }
380  }
381 
382  if ( !found )
383  {
384  return shortName;
385  }
386 
387  int n = 1;
388  shortName = QStringLiteral( "%1_%2" ).arg( fieldName.mid( 0, 8 ) ).arg( n );
389  found = true;
390  while ( found )
391  {
392  found = false;
393  for ( int idx = 0; idx < allFields.count(); ++idx )
394  {
395  if ( shortName == allFields.at( idx ).name() )
396  {
397  n += 1;
398  if ( n < 9 )
399  {
400  shortName = QStringLiteral( "%1_%2" ).arg( fieldName.mid( 0, 8 ) ).arg( n );
401  }
402  else
403  {
404  shortName = QStringLiteral( "%1_%2" ).arg( fieldName.mid( 0, 7 ) ).arg( n );
405  }
406  found = true;
407  }
408  }
409  }
410  return shortName;
411 }
void updateFields()
Will regenerate the fields property of this layer by obtaining all fields from the dataProvider...
QgsFeatureRequest & setDestinationCrs(const QgsCoordinateReferenceSystem &crs, const QgsCoordinateTransformContext &context)
Sets the destination crs for feature&#39;s geometries.
QgsFeatureId id
Definition: qgsfeature.h:64
virtual int bandCount() const =0
Gets number of bands.
Wrapper for iterator of features from vector data provider or vector layer.
A rectangle specified with double values.
Definition: qgsrectangle.h:41
virtual bool addAttributes(const QList< QgsField > &attributes)
Adds new attributes to the provider.
virtual QgsRectangle extent() const
Gets the extent of the interface.
This class provides qgis with the ability to render raster datasets onto the mapcanvas.
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
virtual int ySize() const
QgsWkbTypes::GeometryType geometryType() const
Returns point, line or polygon.
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:122
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
Definition: qgsfeature.h:55
const QgsCoordinateReferenceSystem & crs
bool hasGeometry() const
Returns true if the feature has an associated geometry.
Definition: qgsfeature.cpp:197
Median of pixel values.
Minority of pixel values.
Variety (count of distinct) pixel values.
Base class for feedback objects to be used for cancellation of something running in a worker thread...
Definition: qgsfeedback.h:44
Variance of pixel values.
bool isEmpty() const
Returns true if the rectangle is empty.
Definition: qgsrectangle.h:426
QgsFeatureRequest & setNoAttributes()
Set that no attributes will be fetched.
QgsFields fields() const override=0
Returns the fields associated with this data provider.
virtual bool changeAttributeValues(const QgsChangedAttributesMap &attr_map)
Changes attribute values of existing features.
QgsRectangle intersect(const QgsRectangle &rect) const
Returns the intersection with the given rectangle.
Definition: qgsrectangle.h:312
QMap< int, QVariant > QgsAttributeMap
Definition: qgsattributes.h:38
Sum of pixel values.
Majority of pixel values.
This class wraps a request for features to a vector layer (or directly its vector data provider)...
int fieldNameIndex(const QString &fieldName) const
Returns the index of a field name or -1 if the field does not exist.
Encapsulate a field in an attribute table or data source.
Definition: qgsfield.h:48
int calculateStatistics(QgsFeedback *feedback)
Starts the calculation.
Base class for processing filters like renderers, reprojector, resampler etc.
Max of pixel values.
Mean of pixel values.
Range of pixel values (max - min)
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
QMap< QgsFeatureId, QgsAttributeMap > QgsChangedAttributesMap
Definition: qgsfeature.h:557
QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest()) const override=0
Query the provider for features specified in request.
long featureCount() const override=0
Number of features in the layer.
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.
static QgsProject * instance()
Returns the QgsProject singleton instance.
Definition: qgsproject.cpp:442
QgsRectangle boundingBox() const
Returns the bounding box of the geometry.
This class represents a coordinate reference system (CRS).
QList< QgsField > toList() const
Utility function to return a list of QgsField instances.
Definition: qgsfields.cpp:212
A class that calculates raster statistics (count, sum, mean) for a polygon or multipolygon layer and ...
Min of pixel values.
QgsGeometry geometry
Definition: qgsfeature.h:67
virtual QString storageType() const
Returns the permanent storage type for this layer as a friendly name.
QgsVectorDataProvider * dataProvider() FINAL
Returns the layer&#39;s data provider, it may be nullptr.
bool nextFeature(QgsFeature &f)
This is the base class for vector data providers.
Represents a vector layer which manages a vector based data sets.
virtual int xSize() const
Gets raster size.
Standard deviation of pixel values.