Quantum GIS API Documentation  1.7.5-Wroclaw
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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 
26 #if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800
27 #define TO8(x) (x).toUtf8().constData()
28 #else
29 #define TO8(x) (x).toLocal8Bit().constData()
30 #endif
31 
32 QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer* polygonLayer, const QString& rasterFile, const QString& attributePrefix, int rasterBand )
33  : mRasterFilePath( rasterFile )
34  , mRasterBand( rasterBand )
35  , mPolygonLayer( polygonLayer )
36  , mAttributePrefix( attributePrefix )
37  , mInputNodataValue( -1 )
38 {
39 
40 }
41 
43  : mRasterBand( 0 )
44  , mPolygonLayer( 0 )
45 {
46 
47 }
48 
50 {
51 
52 }
53 
54 int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
55 {
57  {
58  return 1;
59  }
60 
62  if ( !vectorProvider )
63  {
64  return 2;
65  }
66 
67  //open the raster layer and the raster band
68  GDALAllRegister();
69  GDALDatasetH inputDataset = GDALOpen( TO8( mRasterFilePath ), GA_ReadOnly );
70  if ( inputDataset == NULL )
71  {
72  return 3;
73  }
74 
75  if ( GDALGetRasterCount( inputDataset ) < ( mRasterBand - 1 ) )
76  {
77  GDALClose( inputDataset );
78  return 4;
79  }
80 
81  GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset, mRasterBand );
82  if ( rasterBand == NULL )
83  {
84  GDALClose( inputDataset );
85  return 5;
86  }
87  mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, NULL );
88 
89  //get geometry info about raster layer
90  int nCellsX = GDALGetRasterXSize( inputDataset );
91  int nCellsY = GDALGetRasterYSize( inputDataset );
92  double geoTransform[6];
93  if ( GDALGetGeoTransform( inputDataset, geoTransform ) != CE_None )
94  {
95  GDALClose( inputDataset );
96  return 6;
97  }
98  double cellsizeX = geoTransform[1];
99  if ( cellsizeX < 0 )
100  {
101  cellsizeX = -cellsizeX;
102  }
103  double cellsizeY = geoTransform[5];
104  if ( cellsizeY < 0 )
105  {
106  cellsizeY = -cellsizeY;
107  }
108  QgsRectangle rasterBBox( geoTransform[0], geoTransform[3] - ( nCellsY * cellsizeY ), geoTransform[0] + ( nCellsX * cellsizeX ), geoTransform[3] );
109 
110  //add the new count, sum, mean fields to the provider
111  QList<QgsField> newFieldList;
112  QgsField countField( mAttributePrefix + "count", QVariant::Double );
113  QgsField sumField( mAttributePrefix + "sum", QVariant::Double );
114  QgsField meanField( mAttributePrefix + "mean", QVariant::Double );
115  newFieldList.push_back( countField );
116  newFieldList.push_back( sumField );
117  newFieldList.push_back( meanField );
118  vectorProvider->addAttributes( newFieldList );
119 
120  //index of the new fields
121  int countIndex = vectorProvider->fieldNameIndex( mAttributePrefix + "count" );
122  int sumIndex = vectorProvider->fieldNameIndex( mAttributePrefix + "sum" );
123  int meanIndex = vectorProvider->fieldNameIndex( mAttributePrefix + "mean" );
124 
125  if ( countIndex == -1 || sumIndex == -1 || meanIndex == -1 )
126  {
127  return 8;
128  }
129 
130  //progress dialog
131  long featureCount = vectorProvider->featureCount();
132  if ( p )
133  {
134  p->setMaximum( featureCount );
135  }
136 
137 
138  //iterate over each polygon
139  vectorProvider->select( QgsAttributeList(), QgsRectangle(), true, false );
140  QgsFeature f;
141  double count = 0;
142  double sum = 0;
143  double mean = 0;
144  int featureCounter = 0;
145 
146  while ( vectorProvider->nextFeature( f ) )
147  {
148  qWarning( "%d", featureCounter );
149  if ( p )
150  {
151  p->setValue( featureCounter );
152  }
153 
154  if ( p && p->wasCanceled() )
155  {
156  break;
157  }
158 
159  QgsGeometry* featureGeometry = f.geometry();
160  if ( !featureGeometry )
161  {
162  ++featureCounter;
163  continue;
164  }
165 
166  int offsetX, offsetY, nCellsX, nCellsY;
167  if ( cellInfoForBBox( rasterBBox, featureGeometry->boundingBox(), cellsizeX, cellsizeY, offsetX, offsetY, nCellsX, nCellsY ) != 0 )
168  {
169  ++featureCounter;
170  continue;
171  }
172 
173  statisticsFromMiddlePointTest_improved( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
174  rasterBBox, sum, count );
175 
176  if ( count <= 1 )
177  {
178  //the cell resolution is probably larger than the polygon area. We switch to precise pixel - polygon intersection in this case
179  statisticsFromPreciseIntersection( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
180  rasterBBox, sum, count );
181  }
182 
183 
184  if ( count == 0 )
185  {
186  mean = 0;
187  }
188  else
189  {
190  mean = sum / count;
191  }
192 
193  //write the statistics value to the vector data provider
194  QgsChangedAttributesMap changeMap;
195  QgsAttributeMap changeAttributeMap;
196  changeAttributeMap.insert( countIndex, QVariant( count ) );
197  changeAttributeMap.insert( sumIndex, QVariant( sum ) );
198  changeAttributeMap.insert( meanIndex, QVariant( mean ) );
199  changeMap.insert( f.id(), changeAttributeMap );
200  vectorProvider->changeAttributeValues( changeMap );
201 
202  ++featureCounter;
203  }
204 
205  if ( p )
206  {
207  p->setValue( featureCount );
208  }
209 
210  GDALClose( inputDataset );
212 
213  if ( p && p->wasCanceled() )
214  {
215  return 9;
216  }
217 
218  return 0;
219 }
220 
221 int QgsZonalStatistics::cellInfoForBBox( const QgsRectangle& rasterBBox, const QgsRectangle& featureBBox, double cellSizeX, double cellSizeY,
222  int& offsetX, int& offsetY, int& nCellsX, int& nCellsY ) const
223 {
224  //get intersecting bbox
225  QgsRectangle intersectBox = rasterBBox.intersect( &featureBBox );
226  if ( intersectBox.isEmpty() )
227  {
228  nCellsX = 0; nCellsY = 0; offsetX = 0; offsetY = 0;
229  return 0;
230  }
231 
232  //get offset in pixels in x- and y- direction
233  offsetX = ( int )(( intersectBox.xMinimum() - rasterBBox.xMinimum() ) / cellSizeX );
234  offsetY = ( int )(( rasterBBox.yMaximum() - intersectBox.yMaximum() ) / cellSizeY );
235 
236  int maxColumn = ( int )(( intersectBox.xMaximum() - rasterBBox.xMinimum() ) / cellSizeX ) + 1;
237  int maxRow = ( int )(( rasterBBox.yMaximum() - intersectBox.yMinimum() ) / cellSizeY ) + 1;
238 
239  nCellsX = maxColumn - offsetX;
240  nCellsY = maxRow - offsetY;
241 
242  return 0;
243 }
244 
245 void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, QgsGeometry* poly, int pixelOffsetX,
246  int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, double& sum, double& count )
247 {
248  double cellCenterX, cellCenterY;
249  QgsPoint currentCellCenter;
250 
251  float* scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );
252  cellCenterY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
253  count = 0;
254  sum = 0;
255 
256  for ( int i = 0; i < nCellsY; ++i )
257  {
258  GDALRasterIO( band, GF_Read, pixelOffsetX, pixelOffsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 );
259  cellCenterX = rasterBBox.xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;
260  for ( int j = 0; j < nCellsX; ++j )
261  {
262  currentCellCenter = QgsPoint( cellCenterX, cellCenterY );
263  if ( poly->contains( &currentCellCenter ) )
264  {
265  if ( scanLine[j] != mInputNodataValue ) //don't consider nodata values
266  {
267  sum += scanLine[j];
268  ++count;
269  }
270  }
271  cellCenterX += cellSizeX;
272  }
273  cellCenterY -= cellSizeY;
274  }
275  CPLFree( scanLine );
276 }
277 
278 void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, QgsGeometry* poly, int pixelOffsetX,
279  int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, double& sum, double& count )
280 {
281  sum = 0;
282  count = 0;
283  double currentY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
284  float* pixelData = ( float * ) CPLMalloc( sizeof( float ) );
285  QgsGeometry* pixelRectGeometry = 0;
286 
287  double hCellSizeX = cellSizeX / 2.0;
288  double hCellSizeY = cellSizeY / 2.0;
289  double pixelArea = cellSizeX * cellSizeY;
290  double weight = 0;
291 
292  for ( int row = 0; row < nCellsY; ++row )
293  {
294  double currentX = rasterBBox.xMinimum() + cellSizeX / 2.0 + pixelOffsetX * cellSizeX;
295  for ( int col = 0; col < nCellsX; ++col )
296  {
297  GDALRasterIO( band, GF_Read, pixelOffsetX + col, pixelOffsetY + row, nCellsX, 1, pixelData, 1, 1, GDT_Float32, 0, 0 );
298  pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) );
299  if ( pixelRectGeometry )
300  {
301  //intersection
302  QgsGeometry *intersectGeometry = pixelRectGeometry->intersection( poly );
303  if ( intersectGeometry )
304  {
305  double intersectionArea = intersectGeometry->area();
306  if ( intersectionArea >= 0.0 )
307  {
308  weight = intersectionArea / pixelArea;
309  count += weight;
310  sum += *pixelData * weight;
311  }
312  delete intersectGeometry;
313  }
314  }
315  currentX += cellSizeX;
316  }
317  currentY -= cellSizeY;
318  }
319  CPLFree( pixelData );
320 }
321 
322 void QgsZonalStatistics::statisticsFromMiddlePointTest_improved( void* band, QgsGeometry* poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY,
323  double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, double& sum, double& count )
324 {
325  double cellCenterX, cellCenterY;
326  QgsPoint currentCellCenter;
327 
328  float* scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );
329  cellCenterY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
330  count = 0;
331  sum = 0;
332 
333  for ( int i = 0; i < nCellsY; ++i )
334  {
335  GDALRasterIO( band, GF_Read, pixelOffsetX, pixelOffsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 );
336  cellCenterX = rasterBBox.xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;
337 
338  //do intersection of scanline with geometry
339  GEOSCoordSequence* scanLineSequence = GEOSCoordSeq_create( 2, 2 );
340  GEOSCoordSeq_setX( scanLineSequence, 0, cellCenterX );
341  GEOSCoordSeq_setY( scanLineSequence, 0, cellCenterY );
342  GEOSCoordSeq_setX( scanLineSequence, 1, cellCenterX + nCellsX * cellSizeX );
343  GEOSCoordSeq_setY( scanLineSequence, 1, cellCenterY );
344  GEOSGeometry* scanLineGeos = GEOSGeom_createLineString( scanLineSequence ); //todo: delete
345  GEOSGeometry* polyGeos = poly->asGeos();
346  GEOSGeometry* scanLineIntersection = GEOSIntersection( scanLineGeos, polyGeos );
347  GEOSGeom_destroy( scanLineGeos );
348  if ( !scanLineIntersection )
349  {
350  cellCenterY -= cellSizeY;
351  continue;
352  }
353 
354  //debug
355  //char* scanLineIntersectionType = GEOSGeomType( scanLineIntersection );
356 
357  int numGeoms = GEOSGetNumGeometries( scanLineIntersection );
358  if ( numGeoms < 1 )
359  {
360  GEOSGeom_destroy( scanLineIntersection );
361  cellCenterY -= cellSizeY;
362  continue;
363  }
364 
365  QList<double> scanLineList;
366  double currentValue;
367  GEOSGeometry* currentGeom = 0;
368  for ( int z = 0; z < numGeoms; ++z )
369  {
370  if ( numGeoms == 1 )
371  {
372  currentGeom = scanLineIntersection;
373  }
374  else
375  {
376  currentGeom = GEOSGeom_clone( GEOSGetGeometryN( scanLineIntersection, z ) );
377  }
378  const GEOSCoordSequence* scanLineCoordSequence = GEOSGeom_getCoordSeq( currentGeom );
379  if ( !scanLineCoordSequence )
380  {
381  //error
382  }
383  unsigned int scanLineIntersectionSize;
384  GEOSCoordSeq_getSize( scanLineCoordSequence, &scanLineIntersectionSize );
385  if ( !scanLineCoordSequence || scanLineIntersectionSize < 2 || ( scanLineIntersectionSize & 1 ) )
386  {
387  //error
388  }
389  for ( unsigned int k = 0; k < scanLineIntersectionSize; ++k )
390  {
391  GEOSCoordSeq_getX( scanLineCoordSequence, k, &currentValue );
392  scanLineList.push_back( currentValue );
393  }
394 
395  if ( numGeoms != 1 )
396  {
397  GEOSGeom_destroy( currentGeom );
398  }
399  }
400  GEOSGeom_destroy( scanLineIntersection );
401  qSort( scanLineList );
402 
403  if ( scanLineList.size() < 1 )
404  {
405  cellCenterY -= cellSizeY;
406  continue;
407  }
408 
409  int listPlace = -1;
410  for ( int j = 0; j < nCellsX; ++j )
411  {
412  //currentCellCenter = QgsPoint( cellCenterX, cellCenterY );
413 
414  //instead of doing a contained test every time, find the place of scanLineList and check if even / odd
415  if ( listPlace >= scanLineList.size() - 1 )
416  {
417  break;
418  }
419  if ( cellCenterX >= scanLineList.at( listPlace + 1 ) )
420  {
421  ++listPlace;
422  if ( listPlace >= scanLineList.size() )
423  {
424  break;
425  }
426  }
427  if ( listPlace >= 0 && listPlace < ( scanLineList.size() - 1 ) && !( listPlace & 1 ) )
428  {
429  if ( scanLine[j] != mInputNodataValue ) //don't consider nodata values
430  {
431  sum += scanLine[j];
432  ++count;
433  }
434  }
435  cellCenterX += cellSizeX;
436  }
437  cellCenterY -= cellSizeY;
438  }
439  CPLFree( scanLine );
440 }
441 
442