QGIS API Documentation  2.5.0-Master
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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 #include "qgsgeometry.h"
20 #include "qgsvectordataprovider.h"
21 #include "qgsvectorlayer.h"
22 #include "gdal.h"
23 #include "cpl_string.h"
24 #include <QProgressDialog>
25 #include <QFile>
26 
27 #if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800
28 #define TO8F(x) (x).toUtf8().constData()
29 #else
30 #define TO8F(x) QFile::encodeName( x ).constData()
31 #endif
32 
33 QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer* polygonLayer, const QString& rasterFile, const QString& attributePrefix, int rasterBand )
34  : mRasterFilePath( rasterFile )
35  , mRasterBand( rasterBand )
36  , mPolygonLayer( polygonLayer )
37  , mAttributePrefix( attributePrefix )
38  , mInputNodataValue( -1 )
39 {
40 
41 }
42 
44  : mRasterBand( 0 )
45  , mPolygonLayer( 0 )
46 {
47 
48 }
49 
51 {
52 
53 }
54 
55 int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
56 {
58  {
59  return 1;
60  }
61 
63  if ( !vectorProvider )
64  {
65  return 2;
66  }
67 
68  //open the raster layer and the raster band
69  GDALAllRegister();
70  GDALDatasetH inputDataset = GDALOpen( TO8F( mRasterFilePath ), GA_ReadOnly );
71  if ( inputDataset == NULL )
72  {
73  return 3;
74  }
75 
76  if ( GDALGetRasterCount( inputDataset ) < ( mRasterBand - 1 ) )
77  {
78  GDALClose( inputDataset );
79  return 4;
80  }
81 
82  GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset, mRasterBand );
83  if ( rasterBand == NULL )
84  {
85  GDALClose( inputDataset );
86  return 5;
87  }
88  mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, NULL );
89 
90  //get geometry info about raster layer
91  int nCellsXGDAL = GDALGetRasterXSize( inputDataset );
92  int nCellsYGDAL = GDALGetRasterYSize( inputDataset );
93  double geoTransform[6];
94  if ( GDALGetGeoTransform( inputDataset, geoTransform ) != CE_None )
95  {
96  GDALClose( inputDataset );
97  return 6;
98  }
99  double cellsizeX = geoTransform[1];
100  if ( cellsizeX < 0 )
101  {
102  cellsizeX = -cellsizeX;
103  }
104  double cellsizeY = geoTransform[5];
105  if ( cellsizeY < 0 )
106  {
107  cellsizeY = -cellsizeY;
108  }
109  QgsRectangle rasterBBox( geoTransform[0], geoTransform[3] - ( nCellsYGDAL * cellsizeY ),
110  geoTransform[0] + ( nCellsXGDAL * cellsizeX ), geoTransform[3] );
111 
112  //add the new count, sum, mean fields to the provider
113  QList<QgsField> newFieldList;
114  QString countFieldName = getUniqueFieldName( mAttributePrefix + "count" );
115  QString sumFieldName = getUniqueFieldName( mAttributePrefix + "sum" );
116  QString meanFieldName = getUniqueFieldName( mAttributePrefix + "mean" );
117  QgsField countField( countFieldName, QVariant::Double, "double precision" );
118  QgsField sumField( sumFieldName, QVariant::Double, "double precision" );
119  QgsField meanField( meanFieldName, QVariant::Double, "double precision" );
120  newFieldList.push_back( countField );
121  newFieldList.push_back( sumField );
122  newFieldList.push_back( meanField );
123  vectorProvider->addAttributes( newFieldList );
124 
125  //index of the new fields
126  int countIndex = vectorProvider->fieldNameIndex( countFieldName );
127  int sumIndex = vectorProvider->fieldNameIndex( sumFieldName );
128  int meanIndex = vectorProvider->fieldNameIndex( meanFieldName );
129 
130  if ( countIndex == -1 || sumIndex == -1 || meanIndex == -1 )
131  {
132  return 8;
133  }
134 
135  //progress dialog
136  long featureCount = vectorProvider->featureCount();
137  if ( p )
138  {
139  p->setMaximum( featureCount );
140  }
141 
142 
143  //iterate over each polygon
144  QgsFeatureRequest request;
146  QgsFeatureIterator fi = vectorProvider->getFeatures( request );
147  QgsFeature f;
148  double count = 0;
149  double sum = 0;
150  double mean = 0;
151  int featureCounter = 0;
152 
153  while ( fi.nextFeature( f ) )
154  {
155  if ( p )
156  {
157  p->setValue( featureCounter );
158  }
159 
160  if ( p && p->wasCanceled() )
161  {
162  break;
163  }
164 
165  QgsGeometry* featureGeometry = f.geometry();
166  if ( !featureGeometry )
167  {
168  ++featureCounter;
169  continue;
170  }
171 
172  QgsRectangle featureRect = featureGeometry->boundingBox().intersect( &rasterBBox );
173  if ( featureRect.isEmpty() )
174  {
175  ++featureCounter;
176  continue;
177  }
178 
179  int offsetX, offsetY, nCellsX, nCellsY;
180  if ( cellInfoForBBox( rasterBBox, featureRect, cellsizeX, cellsizeY, offsetX, offsetY, nCellsX, nCellsY ) != 0 )
181  {
182  ++featureCounter;
183  continue;
184  }
185 
186  //avoid access to cells outside of the raster (may occur because of rounding)
187  if (( offsetX + nCellsX ) > nCellsXGDAL )
188  {
189  nCellsX = nCellsXGDAL - offsetX;
190  }
191  if (( offsetY + nCellsY ) > nCellsYGDAL )
192  {
193  nCellsY = nCellsYGDAL - offsetY;
194  }
195 
196  statisticsFromMiddlePointTest( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
197  rasterBBox, sum, count );
198 
199  if ( count <= 1 )
200  {
201  //the cell resolution is probably larger than the polygon area. We switch to precise pixel - polygon intersection in this case
202  statisticsFromPreciseIntersection( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
203  rasterBBox, sum, count );
204  }
205 
206 
207  if ( count == 0 )
208  {
209  mean = 0;
210  }
211  else
212  {
213  mean = sum / count;
214  }
215 
216  //write the statistics value to the vector data provider
217  QgsChangedAttributesMap changeMap;
218  QgsAttributeMap changeAttributeMap;
219  changeAttributeMap.insert( countIndex, QVariant( count ) );
220  changeAttributeMap.insert( sumIndex, QVariant( sum ) );
221  changeAttributeMap.insert( meanIndex, QVariant( mean ) );
222  changeMap.insert( f.id(), changeAttributeMap );
223  vectorProvider->changeAttributeValues( changeMap );
224 
225  ++featureCounter;
226  }
227 
228  if ( p )
229  {
230  p->setValue( featureCount );
231  }
232 
233  GDALClose( inputDataset );
235 
236  if ( p && p->wasCanceled() )
237  {
238  return 9;
239  }
240 
241  return 0;
242 }
243 
244 int QgsZonalStatistics::cellInfoForBBox( const QgsRectangle& rasterBBox, const QgsRectangle& featureBBox, double cellSizeX, double cellSizeY,
245  int& offsetX, int& offsetY, int& nCellsX, int& nCellsY ) const
246 {
247  //get intersecting bbox
248  QgsRectangle intersectBox = rasterBBox.intersect( &featureBBox );
249  if ( intersectBox.isEmpty() )
250  {
251  nCellsX = 0; nCellsY = 0; offsetX = 0; offsetY = 0;
252  return 0;
253  }
254 
255  //get offset in pixels in x- and y- direction
256  offsetX = ( int )(( intersectBox.xMinimum() - rasterBBox.xMinimum() ) / cellSizeX );
257  offsetY = ( int )(( rasterBBox.yMaximum() - intersectBox.yMaximum() ) / cellSizeY );
258 
259  int maxColumn = ( int )(( intersectBox.xMaximum() - rasterBBox.xMinimum() ) / cellSizeX ) + 1;
260  int maxRow = ( int )(( rasterBBox.yMaximum() - intersectBox.yMinimum() ) / cellSizeY ) + 1;
261 
262  nCellsX = maxColumn - offsetX;
263  nCellsY = maxRow - offsetY;
264 
265  return 0;
266 }
267 
268 void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, QgsGeometry* poly, int pixelOffsetX,
269  int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, double& sum, double& count )
270 {
271  double cellCenterX, cellCenterY;
272 
273  float* scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );
274  cellCenterY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
275  count = 0;
276  sum = 0;
277 
278  const GEOSGeometry* polyGeos = poly->asGeos();
279  if ( !polyGeos )
280  {
281  return;
282  }
283 
284  const GEOSPreparedGeometry* polyGeosPrepared = GEOSPrepare( poly->asGeos() );
285  if ( !polyGeosPrepared )
286  {
287  return;
288  }
289 
290  GEOSCoordSequence* cellCenterCoords = 0;
291  GEOSGeometry* currentCellCenter = 0;
292 
293  for ( int i = 0; i < nCellsY; ++i )
294  {
295  if ( GDALRasterIO( band, GF_Read, pixelOffsetX, pixelOffsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 )
296  != CPLE_None )
297  {
298  continue;
299  }
300  cellCenterX = rasterBBox.xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;
301  for ( int j = 0; j < nCellsX; ++j )
302  {
303  GEOSGeom_destroy( currentCellCenter );
304  cellCenterCoords = GEOSCoordSeq_create( 1, 2 );
305  GEOSCoordSeq_setX( cellCenterCoords, 0, cellCenterX );
306  GEOSCoordSeq_setY( cellCenterCoords, 0, cellCenterY );
307  currentCellCenter = GEOSGeom_createPoint( cellCenterCoords );
308 
309  if ( scanLine[j] != mInputNodataValue ) //don't consider nodata values
310  {
311  if ( GEOSPreparedContains( polyGeosPrepared, currentCellCenter ) )
312  {
313  if ( !qIsNaN( scanLine[j] ) )
314  {
315  sum += scanLine[j];
316  }
317  ++count;
318  }
319  }
320  cellCenterX += cellSizeX;
321  }
322  cellCenterY -= cellSizeY;
323  }
324  CPLFree( scanLine );
325  GEOSPreparedGeom_destroy( polyGeosPrepared );
326 }
327 
328 void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, QgsGeometry* poly, int pixelOffsetX,
329  int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, double& sum, double& count )
330 {
331  sum = 0;
332  count = 0;
333  double currentY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
334  float* pixelData = ( float * ) CPLMalloc( sizeof( float ) );
335  QgsGeometry* pixelRectGeometry = 0;
336 
337  double hCellSizeX = cellSizeX / 2.0;
338  double hCellSizeY = cellSizeY / 2.0;
339  double pixelArea = cellSizeX * cellSizeY;
340  double weight = 0;
341 
342  for ( int row = 0; row < nCellsY; ++row )
343  {
344  double currentX = rasterBBox.xMinimum() + cellSizeX / 2.0 + pixelOffsetX * cellSizeX;
345  for ( int col = 0; col < nCellsX; ++col )
346  {
347  GDALRasterIO( band, GF_Read, pixelOffsetX + col, pixelOffsetY + row, nCellsX, 1, pixelData, 1, 1, GDT_Float32, 0, 0 );
348  pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) );
349  if ( pixelRectGeometry )
350  {
351  //intersection
352  QgsGeometry *intersectGeometry = pixelRectGeometry->intersection( poly );
353  if ( intersectGeometry )
354  {
355  double intersectionArea = intersectGeometry->area();
356  if ( intersectionArea >= 0.0 )
357  {
358  weight = intersectionArea / pixelArea;
359  count += weight;
360  sum += *pixelData * weight;
361  }
362  delete intersectGeometry;
363  }
364  }
365  currentX += cellSizeX;
366  }
367  currentY -= cellSizeY;
368  }
369  CPLFree( pixelData );
370 }
371 
372 QString QgsZonalStatistics::getUniqueFieldName( QString fieldName )
373 {
375 
376  if ( !dp->storageType().contains( "ESRI Shapefile" ) )
377  {
378  return fieldName;
379  }
380 
381  const QgsFields& providerFields = dp->fields();
382  QString shortName = fieldName.mid( 0, 10 );
383 
384  bool found = false;
385  for ( int idx = 0; idx < providerFields.count(); ++idx )
386  {
387  if ( shortName == providerFields[idx].name() )
388  {
389  found = true;
390  break;
391  }
392  }
393 
394  if ( !found )
395  {
396  return shortName;
397  }
398 
399  int n = 1;
400  shortName = QString( "%1_%2" ).arg( fieldName.mid( 0, 8 ) ).arg( n );
401  found = true;
402  while ( found )
403  {
404  found = false;
405  for ( int idx = 0; idx < providerFields.count(); ++idx )
406  {
407  if ( shortName == providerFields[idx].name() )
408  {
409  n += 1;
410  if ( n < 9 )
411  {
412  shortName = QString( "%1_%2" ).arg( fieldName.mid( 0, 8 ) ).arg( n );
413  }
414  else
415  {
416  shortName = QString( "%1_%2" ).arg( fieldName.mid( 0, 7 ) ).arg( n );
417  }
418  found = true;
419  }
420  }
421  }
422  return shortName;
423 }
QgsFeatureId id() const
Get the feature id for this feature.
Definition: qgsfeature.cpp:100
void updateFields()
Assembles mUpdatedFields considering provider fields, joined fields and added fields.
Wrapper for iterator of features from vector data provider or vector layer.
A rectangle specified with double values.
Definition: qgsrectangle.h:35
bool isEmpty() const
test if rectangle is empty.
QMap< int, QVariant > QgsAttributeMap
Definition: qgsfeature.h:98
virtual bool addAttributes(const QList< QgsField > &attributes)
Adds new attributes.
double yMaximum() const
Get the y maximum value (top side of rectangle)
Definition: qgsrectangle.h:194
QgsGeometry * geometry() const
Get the geometry object associated with this feature.
Definition: qgsfeature.cpp:112
QgsFeatureRequest & setSubsetOfAttributes(const QgsAttributeList &attrs)
Set a subset of attributes that will be fetched.
Container of fields for a vector layer.
Definition: qgsfield.h:163
void statisticsFromMiddlePointTest(void *band, QgsGeometry *poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle &rasterBBox, double &sum, double &count)
Returns statistics by considering the pixels where the center point is within the polygon (fast) ...
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
Definition: qgsfeature.h:113
int fieldNameIndex(const QString &fieldName) const
Returns the index of a field name or -1 if the field does not exist.
double area()
get area of geometry using GEOS
void statisticsFromPreciseIntersection(void *band, QgsGeometry *poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle &rasterBBox, double &sum, double &count)
Returns statistics with precise pixel - polygon intersection test (slow)
QString getUniqueFieldName(QString fieldName)
const GEOSGeometry * asGeos() const
Returns a geos geomtry.
int mRasterBand
Raster band to calculate statistics from (defaults to 1)
double yMinimum() const
Get the y minimum value (bottom side of rectangle)
Definition: qgsrectangle.h:199
double xMaximum() const
Get the x maximum value (right side of rectangle)
Definition: qgsrectangle.h:184
virtual bool changeAttributeValues(const QgsChangedAttributesMap &attr_map)
Changes attribute values of existing features.
virtual QString storageType() const
Returns the permanent storage type for this layer as a friendly name.
virtual long featureCount() const =0
Number of features in the layer.
This class wraps a request for features to a vector layer (or directly its vector data provider)...
QList< int > QgsAttributeList
virtual QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest())=0
Query the provider for features specified in request.
int count() const
Return number of items.
Definition: qgsfield.h:200
QGis::GeometryType geometryType() const
Returns point, line or polygon.
Encapsulate a field in an attribute table or data source.
Definition: qgsfield.h:33
QgsGeometry * intersection(QgsGeometry *geometry)
Returns a geometry representing the points shared by this geometry and other.
int cellInfoForBBox(const QgsRectangle &rasterBBox, const QgsRectangle &featureBBox, double cellSizeX, double cellSizeY, int &offsetX, int &offsetY, int &nCellsX, int &nCellsY) const
Analysis what cells need to be considered to cover the bounding box of a feature. ...
#define TO8F(x)
QgsRectangle boundingBox()
Returns the bounding box of this feature.
virtual const QgsFields & fields() const =0
Return a map of indexes with field names for this layer.
int calculateStatistics(QProgressDialog *p)
Starts the calculation.
QMap< QgsFeatureId, QgsAttributeMap > QgsChangedAttributesMap
Definition: qgsfeature.h:320
QgsVectorLayer * mPolygonLayer
QgsRectangle intersect(const QgsRectangle *rect) const
return the intersection with the given rectangle
static QgsGeometry * fromRect(const QgsRectangle &rect)
construct geometry from a rectangle
float mInputNodataValue
The nodata value of the input layer.
QgsVectorDataProvider * dataProvider()
Returns the data provider.
bool nextFeature(QgsFeature &f)
This is the base class for vector data providers.
Represents a vector layer which manages a vector based data sets.
double xMinimum() const
Get the x minimum value (left side of rectangle)
Definition: qgsrectangle.h:189