← Back to branch summary

~ubuntu-branches/ubuntu/trusty/qgis/trusty

« back to all changes in this revision

Viewing changes to src/analysis/vector/qgszonalstatistics.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Johan Van de Wauw
  • Date: 2010-07-11 20:23:24 UTC
  • mfrom: (3.1.4 squeeze)
  • Revision ID: james.westby@ubuntu.com-20100711202324-5ktghxa7hracohmr
Tags: 1.4.0+12730-3ubuntu1
https://launchpad.net/bugs/540941
* Merge from Debian unstable (LP: #540941).
* Fix compilation issues with QT 4.7
* Add build-depends on libqt4-webkit-dev 

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/analysis/vector/qgszonalstatistics.cpp
 
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
QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer* polygonLayer, const QString& rasterFile, const QString& attributePrefix, int rasterBand ) :
 
27
    mRasterFilePath( rasterFile ), mRasterBand( rasterBand ), mPolygonLayer( polygonLayer ), mAttributePrefix( attributePrefix ), mInputNodataValue( -1 )
 
28
{
 
29
 
 
30
}
 
31
 
 
32
QgsZonalStatistics::QgsZonalStatistics(): mRasterBand( 0 ), mPolygonLayer( 0 )
 
33
{
 
34
 
 
35
}
 
36
 
 
37
QgsZonalStatistics::~QgsZonalStatistics()
 
38
{
 
39
 
 
40
}
 
41
 
 
42
int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
 
43
{
 
44
  if ( !mPolygonLayer || mPolygonLayer->geometryType() != QGis::Polygon )
 
45
  {
 
46
    return 1;
 
47
  }
 
48
 
 
49
  QgsVectorDataProvider* vectorProvider = mPolygonLayer->dataProvider();
 
50
  if ( !vectorProvider )
 
51
  {
 
52
    return 2;
 
53
  }
 
54
 
 
55
  //open the raster layer and the raster band
 
56
  GDALAllRegister();
 
57
  GDALDatasetH inputDataset = GDALOpen( mRasterFilePath.toLocal8Bit().data(), GA_ReadOnly );
 
58
  if ( inputDataset == NULL )
 
59
  {
 
60
    return 3;
 
61
  }
 
62
 
 
63
  if ( GDALGetRasterCount( inputDataset ) < ( mRasterBand - 1 ) )
 
64
  {
 
65
    GDALClose( inputDataset );
 
66
    return 4;
 
67
  }
 
68
 
 
69
  GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset, mRasterBand );
 
70
  if ( rasterBand == NULL )
 
71
  {
 
72
    GDALClose( inputDataset );
 
73
    return 5;
 
74
  }
 
75
  mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, NULL );
 
76
 
 
77
  //get geometry info about raster layer
 
78
  int nCellsX = GDALGetRasterXSize( inputDataset );
 
79
  int nCellsY = GDALGetRasterYSize( inputDataset );
 
80
  double geoTransform[6];
 
81
  if ( GDALGetGeoTransform( inputDataset, geoTransform ) != CE_None )
 
82
  {
 
83
    GDALClose( inputDataset );
 
84
    return 6;
 
85
  }
 
86
  double cellsizeX = geoTransform[1];
 
87
  if ( cellsizeX < 0 )
 
88
  {
 
89
    cellsizeX = -cellsizeX;
 
90
  }
 
91
  double cellsizeY = geoTransform[5];
 
92
  if ( cellsizeY < 0 )
 
93
  {
 
94
    cellsizeY = -cellsizeY;
 
95
  }
 
96
  QgsRectangle rasterBBox( geoTransform[0], geoTransform[3] - ( nCellsY * cellsizeY ), geoTransform[0] + ( nCellsX * cellsizeX ), geoTransform[3] );
 
97
 
 
98
  //add the new count, sum, mean fields to the provider
 
99
  QList<QgsField> newFieldList;
 
100
  QgsField countField( mAttributePrefix + "count", QVariant::Int );
 
101
  QgsField sumField( mAttributePrefix + "sum", QVariant::Double );
 
102
  QgsField meanField( mAttributePrefix + "mean", QVariant::Double );
 
103
  newFieldList.push_back( countField );
 
104
  newFieldList.push_back( sumField );
 
105
  newFieldList.push_back( meanField );
 
106
  if ( !vectorProvider->addAttributes( newFieldList ) )
 
107
  {
 
108
    return 7;
 
109
  }
 
110
 
 
111
  //index of the new fields
 
112
  int countIndex = vectorProvider->fieldNameIndex( mAttributePrefix + "count" );
 
113
  int sumIndex = vectorProvider->fieldNameIndex( mAttributePrefix + "sum" );
 
114
  int meanIndex = vectorProvider->fieldNameIndex( mAttributePrefix + "mean" );
 
115
 
 
116
  if ( countIndex == -1 || sumIndex == -1 || meanIndex == -1 )
 
117
  {
 
118
    return 8;
 
119
  }
 
120
 
 
121
  //progress dialog
 
122
  long featureCount = vectorProvider->featureCount();
 
123
  if ( p )
 
124
  {
 
125
    p->setMaximum( featureCount );
 
126
  }
 
127
 
 
128
 
 
129
  //iterate over each polygon
 
130
  vectorProvider->select( QgsAttributeList(), QgsRectangle(), true, false );
 
131
  vectorProvider->rewind();
 
132
  QgsFeature f;
 
133
  double count = 0;
 
134
  double sum = 0;
 
135
  double mean = 0;
 
136
  float* scanLine;
 
137
  int featureCounter = 0;
 
138
  //x- and y- coordinate of current cell
 
139
  double cellCenterX = 0;
 
140
  double cellCenterY = 0;
 
141
  QgsPoint currentCellCenter;
 
142
 
 
143
  while ( vectorProvider->nextFeature( f ) )
 
144
  {
 
145
    if ( p )
 
146
    {
 
147
      p->setValue( featureCounter );
 
148
    }
 
149
 
 
150
    if ( p && p->wasCanceled() )
 
151
    {
 
152
      break;
 
153
    }
 
154
 
 
155
    QgsGeometry* featureGeometry = f.geometry();
 
156
    if ( !featureGeometry )
 
157
    {
 
158
      ++featureCounter;
 
159
      continue;
 
160
    }
 
161
 
 
162
    int offsetX, offsetY, nCellsX, nCellsY;
 
163
    if ( cellInfoForBBox( rasterBBox, featureGeometry->boundingBox(), cellsizeX, cellsizeY, offsetX, offsetY, nCellsX, nCellsY ) != 0 )
 
164
    {
 
165
      ++featureCounter;
 
166
      continue;
 
167
    }
 
168
 
 
169
    scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );
 
170
    cellCenterY = rasterBBox.yMaximum() - offsetY * cellsizeY - cellsizeY / 2;
 
171
    count = 0;
 
172
    sum = 0;
 
173
 
 
174
    for ( int i = 0; i < nCellsY; ++i )
 
175
    {
 
176
      GDALRasterIO( rasterBand, GF_Read, offsetX, offsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 );
 
177
      cellCenterX = rasterBBox.xMinimum() + offsetX * cellsizeX + cellsizeX / 2;
 
178
      for ( int j = 0; j < nCellsX; ++j )
 
179
      {
 
180
        currentCellCenter = QgsPoint( cellCenterX, cellCenterY );
 
181
        if ( featureGeometry->contains( &currentCellCenter ) )
 
182
        {
 
183
          if ( scanLine[j] != mInputNodataValue ) //don't consider nodata values
 
184
          {
 
185
            sum += scanLine[j];
 
186
            ++count;
 
187
          }
 
188
        }
 
189
        cellCenterX += cellsizeX;
 
190
      }
 
191
      cellCenterY -= cellsizeY;
 
192
    }
 
193
 
 
194
    if ( count == 0 )
 
195
    {
 
196
      mean = 0;
 
197
    }
 
198
    else
 
199
    {
 
200
      mean = sum / count;
 
201
    }
 
202
 
 
203
    //write the new AEY value to the vector data provider
 
204
    QgsChangedAttributesMap changeMap;
 
205
    QgsAttributeMap changeAttributeMap;
 
206
    changeAttributeMap.insert( countIndex, QVariant( count ) );
 
207
    changeAttributeMap.insert( sumIndex, QVariant( sum ) );
 
208
    changeAttributeMap.insert( meanIndex, QVariant( mean ) );
 
209
    changeMap.insert( f.id(), changeAttributeMap );
 
210
    vectorProvider->changeAttributeValues( changeMap );
 
211
 
 
212
    CPLFree( scanLine );
 
213
    ++featureCounter;
 
214
  }
 
215
 
 
216
  if ( p )
 
217
  {
 
218
    p->setValue( featureCount );
 
219
  }
 
220
 
 
221
  GDALClose( inputDataset );
 
222
  return 0;
 
223
}
 
224
 
 
225
int QgsZonalStatistics::cellInfoForBBox( const QgsRectangle& rasterBBox, const QgsRectangle& featureBBox, double cellSizeX, double cellSizeY,
 
226
    int& offsetX, int& offsetY, int& nCellsX, int& nCellsY ) const
 
227
{
 
228
  //get intersecting bbox
 
229
  QgsRectangle intersectBox = rasterBBox.intersect( &featureBBox );
 
230
  if ( intersectBox.isEmpty() )
 
231
  {
 
232
    nCellsX = 0; nCellsY = 0; offsetX = 0; offsetY = 0;
 
233
    return 0;
 
234
  }
 
235
 
 
236
  //get offset in pixels in x- and y- direction
 
237
  offsetX = ( int )(( intersectBox.xMinimum() - rasterBBox.xMinimum() ) / cellSizeX );
 
238
  offsetY = ( int )(( rasterBBox.yMaximum() - intersectBox.yMaximum() ) / cellSizeY );
 
239
 
 
240
  int maxColumn = ( int )(( intersectBox.xMaximum() - rasterBBox.xMinimum() ) / cellSizeX ) + 1;
 
241
  int maxRow = ( int )(( rasterBBox.yMaximum() - intersectBox.yMinimum() ) / cellSizeY ) + 1;
 
242
 
 
243
  nCellsX = maxColumn - offsetX;
 
244
  nCellsY = maxRow - offsetY;
 
245
 
 
246
  return 0;
 
247
}
 
248
 
 
249
 
 
250
 
 
251