Branch data Line data Source code
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"
25 : : #include "processing/qgsrasteranalysisutils.h"
26 : : #include "qgsrasterdataprovider.h"
27 : : #include "qgsrasterlayer.h"
28 : : #include "qgslogger.h"
29 : : #include "qgsproject.h"
30 : :
31 : : #include <QFile>
32 : :
33 : 0 : QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer *polygonLayer, QgsRasterLayer *rasterLayer, const QString &attributePrefix, int rasterBand, QgsZonalStatistics::Statistics stats )
34 : 0 : : QgsZonalStatistics( polygonLayer,
35 : 0 : rasterLayer ? rasterLayer->dataProvider() : nullptr,
36 : 0 : rasterLayer ? rasterLayer->crs() : QgsCoordinateReferenceSystem(),
37 : 0 : rasterLayer ? rasterLayer->rasterUnitsPerPixelX() : 0,
38 : 0 : rasterLayer ? rasterLayer->rasterUnitsPerPixelY() : 0,
39 : 0 : attributePrefix,
40 : 0 : rasterBand,
41 : 0 : stats )
42 : : {
43 : 0 : }
44 : :
45 : 0 : QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer *polygonLayer, QgsRasterInterface *rasterInterface,
46 : : const QgsCoordinateReferenceSystem &rasterCrs, double rasterUnitsPerPixelX, double rasterUnitsPerPixelY, const QString &attributePrefix, int rasterBand, QgsZonalStatistics::Statistics stats )
47 : 0 : : mRasterInterface( rasterInterface )
48 : 0 : , mRasterCrs( rasterCrs )
49 : 0 : , mCellSizeX( std::fabs( rasterUnitsPerPixelX ) )
50 : 0 : , mCellSizeY( std::fabs( rasterUnitsPerPixelY ) )
51 : 0 : , mRasterBand( rasterBand )
52 : 0 : , mPolygonLayer( polygonLayer )
53 : 0 : , mAttributePrefix( attributePrefix )
54 : 0 : , mStatistics( stats )
55 : : {
56 : 0 : }
57 : :
58 : 0 : QgsZonalStatistics::Result QgsZonalStatistics::calculateStatistics( QgsFeedback *feedback )
59 : : {
60 : 0 : if ( !mRasterInterface )
61 : : {
62 : 0 : return RasterInvalid;
63 : : }
64 : :
65 : 0 : if ( mRasterInterface->bandCount() < mRasterBand )
66 : : {
67 : 0 : return RasterBandInvalid;
68 : : }
69 : :
70 : 0 : if ( !mPolygonLayer || mPolygonLayer->geometryType() != QgsWkbTypes::PolygonGeometry )
71 : : {
72 : 0 : return LayerTypeWrong;
73 : : }
74 : :
75 : 0 : QgsVectorDataProvider *vectorProvider = mPolygonLayer->dataProvider();
76 : 0 : if ( !vectorProvider )
77 : : {
78 : 0 : return LayerInvalid;
79 : : }
80 : :
81 : 0 : QMap<QgsZonalStatistics::Statistic, int> statFieldIndexes;
82 : :
83 : : //add the new fields to the provider
84 : 0 : int oldFieldCount = vectorProvider->fields().count();
85 : 0 : QList<QgsField> newFieldList;
86 : 0 : for ( QgsZonalStatistics::Statistic stat :
87 : 0 : {
88 : : QgsZonalStatistics::Count,
89 : : QgsZonalStatistics::Sum,
90 : : QgsZonalStatistics::Mean,
91 : : QgsZonalStatistics::Median,
92 : : QgsZonalStatistics::StDev,
93 : : QgsZonalStatistics::Min,
94 : : QgsZonalStatistics::Max,
95 : : QgsZonalStatistics::Range,
96 : : QgsZonalStatistics::Minority,
97 : : QgsZonalStatistics::Majority,
98 : : QgsZonalStatistics::Variety,
99 : : QgsZonalStatistics::Variance
100 : : } )
101 : : {
102 : 0 : if ( mStatistics & stat )
103 : : {
104 : 0 : QString fieldName = getUniqueFieldName( mAttributePrefix + QgsZonalStatistics::shortName( stat ), newFieldList );
105 : 0 : QgsField field( fieldName, QVariant::Double, QStringLiteral( "double precision" ) );
106 : 0 : newFieldList.push_back( field );
107 : 0 : statFieldIndexes.insert( stat, oldFieldCount + newFieldList.count() - 1 );
108 : 0 : }
109 : : }
110 : :
111 : 0 : vectorProvider->addAttributes( newFieldList );
112 : :
113 : 0 : long featureCount = vectorProvider->featureCount();
114 : :
115 : 0 : QgsFeatureRequest request;
116 : 0 : request.setNoAttributes();
117 : :
118 : 0 : request.setDestinationCrs( mRasterCrs, QgsProject::instance()->transformContext() );
119 : 0 : QgsFeatureIterator fi = vectorProvider->getFeatures( request );
120 : 0 : QgsFeature feature;
121 : :
122 : 0 : int featureCounter = 0;
123 : :
124 : 0 : QgsChangedAttributesMap changeMap;
125 : 0 : while ( fi.nextFeature( feature ) )
126 : : {
127 : 0 : ++featureCounter;
128 : 0 : if ( feedback && feedback->isCanceled() )
129 : : {
130 : 0 : break;
131 : : }
132 : :
133 : 0 : if ( feedback )
134 : : {
135 : 0 : feedback->setProgress( 100.0 * static_cast< double >( featureCounter ) / featureCount );
136 : 0 : }
137 : :
138 : 0 : QgsGeometry featureGeometry = feature.geometry();
139 : :
140 : 0 : QMap<QgsZonalStatistics::Statistic, QVariant> results = calculateStatistics( mRasterInterface, featureGeometry, mCellSizeX, mCellSizeY, mRasterBand, mStatistics );
141 : :
142 : 0 : if ( results.empty() )
143 : 0 : continue;
144 : :
145 : 0 : QgsAttributeMap changeAttributeMap;
146 : 0 : for ( const auto &result : results.toStdMap() )
147 : : {
148 : 0 : changeAttributeMap.insert( statFieldIndexes.value( result.first ), result.second );
149 : : }
150 : :
151 : 0 : changeMap.insert( feature.id(), changeAttributeMap );
152 : 0 : }
153 : :
154 : 0 : vectorProvider->changeAttributeValues( changeMap );
155 : 0 : mPolygonLayer->updateFields();
156 : :
157 : 0 : if ( feedback )
158 : : {
159 : 0 : if ( feedback->isCanceled() )
160 : 0 : return Canceled;
161 : :
162 : 0 : feedback->setProgress( 100 );
163 : 0 : }
164 : :
165 : 0 : return Success;
166 : 0 : }
167 : :
168 : 0 : QString QgsZonalStatistics::getUniqueFieldName( const QString &fieldName, const QList<QgsField> &newFields )
169 : : {
170 : 0 : QgsVectorDataProvider *dp = mPolygonLayer->dataProvider();
171 : :
172 : 0 : if ( !dp->storageType().contains( QLatin1String( "ESRI Shapefile" ) ) )
173 : : {
174 : 0 : return fieldName;
175 : : }
176 : :
177 : 0 : QList<QgsField> allFields = dp->fields().toList();
178 : 0 : allFields.append( newFields );
179 : 0 : QString shortName = fieldName.mid( 0, 10 );
180 : :
181 : 0 : bool found = false;
182 : 0 : for ( const QgsField &field : std::as_const( allFields ) )
183 : : {
184 : 0 : if ( shortName == field.name() )
185 : : {
186 : 0 : found = true;
187 : 0 : break;
188 : : }
189 : : }
190 : :
191 : 0 : if ( !found )
192 : : {
193 : 0 : return shortName;
194 : : }
195 : :
196 : 0 : int n = 1;
197 : 0 : shortName = QStringLiteral( "%1_%2" ).arg( fieldName.mid( 0, 8 ) ).arg( n );
198 : 0 : found = true;
199 : 0 : while ( found )
200 : : {
201 : 0 : found = false;
202 : 0 : for ( const QgsField &field : std::as_const( allFields ) )
203 : : {
204 : 0 : if ( shortName == field.name() )
205 : : {
206 : 0 : n += 1;
207 : 0 : if ( n < 9 )
208 : : {
209 : 0 : shortName = QStringLiteral( "%1_%2" ).arg( fieldName.mid( 0, 8 ) ).arg( n );
210 : 0 : }
211 : : else
212 : : {
213 : 0 : shortName = QStringLiteral( "%1_%2" ).arg( fieldName.mid( 0, 7 ) ).arg( n );
214 : : }
215 : 0 : found = true;
216 : 0 : }
217 : : }
218 : : }
219 : 0 : return shortName;
220 : 0 : }
221 : :
222 : 0 : QString QgsZonalStatistics::displayName( QgsZonalStatistics::Statistic statistic )
223 : : {
224 : 0 : switch ( statistic )
225 : : {
226 : : case Count:
227 : 0 : return QObject::tr( "Count" );
228 : : case Sum:
229 : 0 : return QObject::tr( "Sum" );
230 : : case Mean:
231 : 0 : return QObject::tr( "Mean" );
232 : : case Median:
233 : 0 : return QObject::tr( "Median" );
234 : : case StDev:
235 : 0 : return QObject::tr( "St dev" );
236 : : case Min:
237 : 0 : return QObject::tr( "Minimum" );
238 : : case Max:
239 : 0 : return QObject::tr( "Maximum" );
240 : : case Range:
241 : 0 : return QObject::tr( "Range" );
242 : : case Minority:
243 : 0 : return QObject::tr( "Minority" );
244 : : case Majority:
245 : 0 : return QObject::tr( "Majority" );
246 : : case Variety:
247 : 0 : return QObject::tr( "Variety" );
248 : : case Variance:
249 : 0 : return QObject::tr( "Variance" );
250 : : case All:
251 : 0 : return QString();
252 : : }
253 : 0 : return QString();
254 : 0 : }
255 : :
256 : 0 : QString QgsZonalStatistics::shortName( QgsZonalStatistics::Statistic statistic )
257 : : {
258 : 0 : switch ( statistic )
259 : : {
260 : : case Count:
261 : 0 : return QStringLiteral( "count" );
262 : : case Sum:
263 : 0 : return QStringLiteral( "sum" );
264 : : case Mean:
265 : 0 : return QStringLiteral( "mean" );
266 : : case Median:
267 : 0 : return QStringLiteral( "median" );
268 : : case StDev:
269 : 0 : return QStringLiteral( "stdev" );
270 : : case Min:
271 : 0 : return QStringLiteral( "min" );
272 : : case Max:
273 : 0 : return QStringLiteral( "max" );
274 : : case Range:
275 : 0 : return QStringLiteral( "range" );
276 : : case Minority:
277 : 0 : return QStringLiteral( "minority" );
278 : : case Majority:
279 : 0 : return QStringLiteral( "majority" );
280 : : case Variety:
281 : 0 : return QStringLiteral( "variety" );
282 : : case Variance:
283 : 0 : return QStringLiteral( "variance" );
284 : : case All:
285 : 0 : return QString();
286 : : }
287 : 0 : return QString();
288 : 0 : }
289 : :
290 : 0 : QMap<QgsZonalStatistics::Statistic, QVariant> QgsZonalStatistics::calculateStatistics( QgsRasterInterface *rasterInterface, const QgsGeometry &geometry, double cellSizeX, double cellSizeY, int rasterBand, QgsZonalStatistics::Statistics statistics )
291 : : {
292 : 0 : QMap<QgsZonalStatistics::Statistic, QVariant> results;
293 : :
294 : 0 : if ( geometry.isEmpty() )
295 : 0 : return results;
296 : :
297 : 0 : QgsRectangle rasterBBox = rasterInterface->extent();
298 : :
299 : 0 : QgsRectangle featureRect = geometry.boundingBox().intersect( rasterBBox );
300 : :
301 : 0 : if ( featureRect.isEmpty() )
302 : 0 : return results;
303 : :
304 : 0 : bool statsStoreValues = ( statistics & QgsZonalStatistics::Median ) ||
305 : 0 : ( statistics & QgsZonalStatistics::StDev ) ||
306 : 0 : ( statistics & QgsZonalStatistics::Variance );
307 : 0 : bool statsStoreValueCount = ( statistics & QgsZonalStatistics::Minority ) ||
308 : 0 : ( statistics & QgsZonalStatistics::Majority );
309 : :
310 : 0 : FeatureStats featureStats( statsStoreValues, statsStoreValueCount );
311 : :
312 : 0 : int nCellsXProvider = rasterInterface->xSize();
313 : 0 : int nCellsYProvider = rasterInterface->ySize();
314 : :
315 : : int nCellsX, nCellsY;
316 : 0 : QgsRectangle rasterBlockExtent;
317 : 0 : QgsRasterAnalysisUtils::cellInfoForBBox( rasterBBox, featureRect, cellSizeX, cellSizeY, nCellsX, nCellsY, nCellsXProvider, nCellsYProvider, rasterBlockExtent );
318 : :
319 : 0 : featureStats.reset();
320 : 0 : QgsRasterAnalysisUtils::statisticsFromMiddlePointTest( rasterInterface, rasterBand, geometry, nCellsX, nCellsY, cellSizeX, cellSizeY, rasterBlockExtent, [ &featureStats ]( double value ) { featureStats.addValue( value ); } );
321 : :
322 : 0 : if ( featureStats.count <= 1 )
323 : : {
324 : : //the cell resolution is probably larger than the polygon area. We switch to precise pixel - polygon intersection in this case
325 : 0 : featureStats.reset();
326 : 0 : QgsRasterAnalysisUtils::statisticsFromPreciseIntersection( rasterInterface, rasterBand, geometry, nCellsX, nCellsY, cellSizeX, cellSizeY, rasterBlockExtent, [ &featureStats ]( double value, double weight ) { featureStats.addValue( value, weight ); } );
327 : 0 : }
328 : :
329 : : // calculate the statistics
330 : 0 : QgsAttributeMap changeAttributeMap;
331 : 0 : if ( statistics & QgsZonalStatistics::Count )
332 : 0 : results.insert( QgsZonalStatistics::Count, QVariant( featureStats.count ) );
333 : 0 : if ( statistics & QgsZonalStatistics::Sum )
334 : 0 : results.insert( QgsZonalStatistics::Sum, QVariant( featureStats.sum ) );
335 : 0 : if ( featureStats.count > 0 )
336 : : {
337 : 0 : double mean = featureStats.sum / featureStats.count;
338 : 0 : if ( statistics & QgsZonalStatistics::Mean )
339 : 0 : results.insert( QgsZonalStatistics::Mean, QVariant( mean ) );
340 : 0 : if ( statistics & QgsZonalStatistics::Median )
341 : : {
342 : 0 : std::sort( featureStats.values.begin(), featureStats.values.end() );
343 : 0 : int size = featureStats.values.count();
344 : 0 : bool even = ( size % 2 ) < 1;
345 : : double medianValue;
346 : 0 : if ( even )
347 : : {
348 : 0 : medianValue = ( featureStats.values.at( size / 2 - 1 ) + featureStats.values.at( size / 2 ) ) / 2;
349 : 0 : }
350 : : else //odd
351 : : {
352 : 0 : medianValue = featureStats.values.at( ( size + 1 ) / 2 - 1 );
353 : : }
354 : 0 : results.insert( QgsZonalStatistics::Median, QVariant( medianValue ) );
355 : 0 : }
356 : 0 : if ( statistics & QgsZonalStatistics::StDev || statistics & QgsZonalStatistics::Variance )
357 : : {
358 : 0 : double sumSquared = 0;
359 : 0 : for ( int i = 0; i < featureStats.values.count(); ++i )
360 : : {
361 : 0 : double diff = featureStats.values.at( i ) - mean;
362 : 0 : sumSquared += diff * diff;
363 : 0 : }
364 : 0 : double variance = sumSquared / featureStats.values.count();
365 : 0 : if ( statistics & QgsZonalStatistics::StDev )
366 : : {
367 : 0 : double stdev = std::pow( variance, 0.5 );
368 : 0 : results.insert( QgsZonalStatistics::StDev, QVariant( stdev ) );
369 : 0 : }
370 : 0 : if ( statistics & QgsZonalStatistics::Variance )
371 : 0 : results.insert( QgsZonalStatistics::Variance, QVariant( variance ) );
372 : 0 : }
373 : 0 : if ( statistics & QgsZonalStatistics::Min )
374 : 0 : results.insert( QgsZonalStatistics::Min, QVariant( featureStats.min ) );
375 : 0 : if ( statistics & QgsZonalStatistics::Max )
376 : 0 : results.insert( QgsZonalStatistics::Max, QVariant( featureStats.max ) );
377 : 0 : if ( statistics & QgsZonalStatistics::Range )
378 : 0 : results.insert( QgsZonalStatistics::Range, QVariant( featureStats.max - featureStats.min ) );
379 : 0 : if ( statistics & QgsZonalStatistics::Minority || statistics & QgsZonalStatistics::Majority )
380 : : {
381 : 0 : QList<int> vals = featureStats.valueCount.values();
382 : 0 : std::sort( vals.begin(), vals.end() );
383 : 0 : if ( statistics & QgsZonalStatistics::Minority )
384 : : {
385 : 0 : double minorityKey = featureStats.valueCount.key( vals.first() );
386 : 0 : results.insert( QgsZonalStatistics::Minority, QVariant( minorityKey ) );
387 : 0 : }
388 : 0 : if ( statistics & QgsZonalStatistics::Majority )
389 : : {
390 : 0 : double majKey = featureStats.valueCount.key( vals.last() );
391 : 0 : results.insert( QgsZonalStatistics::Majority, QVariant( majKey ) );
392 : 0 : }
393 : 0 : }
394 : 0 : if ( statistics & QgsZonalStatistics::Variety )
395 : 0 : results.insert( QgsZonalStatistics::Variety, QVariant( featureStats.valueCount.count() ) );
396 : 0 : }
397 : :
398 : 0 : return results;
399 : 0 : }
|
