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