24 #include "cpl_string.h"
25 #include <QProgressDialog>
28 #if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800
29 #define TO8F(x) (x).toUtf8().constData()
31 #define TO8F(x) QFile::encodeName( x ).constData()
35 : mRasterFilePath( rasterFile )
36 , mRasterBand( rasterBand )
37 , mPolygonLayer( polygonLayer )
38 , mAttributePrefix( attributePrefix )
39 , mInputNodataValue( -1 )
40 , mStatistics( stats )
45 QgsZonalStatistics::QgsZonalStatistics()
48 , mInputNodataValue( -1 )
67 if ( !vectorProvider )
74 GDALDatasetH inputDataset = GDALOpen(
TO8F( mRasterFilePath ), GA_ReadOnly );
75 if ( inputDataset == NULL )
80 if ( GDALGetRasterCount( inputDataset ) < ( mRasterBand - 1 ) )
82 GDALClose( inputDataset );
86 GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset, mRasterBand );
87 if ( rasterBand == NULL )
89 GDALClose( inputDataset );
92 mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, NULL );
95 int nCellsXGDAL = GDALGetRasterXSize( inputDataset );
96 int nCellsYGDAL = GDALGetRasterYSize( inputDataset );
97 double geoTransform[6];
98 if ( GDALGetGeoTransform( inputDataset, geoTransform ) != CE_None )
100 GDALClose( inputDataset );
103 double cellsizeX = geoTransform[1];
106 cellsizeX = -cellsizeX;
108 double cellsizeY = geoTransform[5];
111 cellsizeY = -cellsizeY;
113 QgsRectangle rasterBBox( geoTransform[0], geoTransform[3] - ( nCellsYGDAL * cellsizeY ),
114 geoTransform[0] + ( nCellsXGDAL * cellsizeX ), geoTransform[3] );
121 countFieldName = getUniqueFieldName( mAttributePrefix +
"count" );
122 QgsField countField( countFieldName, QVariant::Double,
"double precision" );
128 sumFieldName = getUniqueFieldName( mAttributePrefix +
"sum" );
129 QgsField sumField( sumFieldName, QVariant::Double,
"double precision" );
135 meanFieldName = getUniqueFieldName( mAttributePrefix +
"mean" );
136 QgsField meanField( meanFieldName, QVariant::Double,
"double precision" );
142 medianFieldName = getUniqueFieldName( mAttributePrefix +
"median" );
143 QgsField medianField( medianFieldName, QVariant::Double,
"double precision" );
149 stdevFieldName = getUniqueFieldName( mAttributePrefix +
"stdev" );
150 QgsField stdField( stdevFieldName, QVariant::Double,
"double precision" );
156 minFieldName = getUniqueFieldName( mAttributePrefix +
"min" );
157 QgsField minField( minFieldName, QVariant::Double,
"double precision" );
163 maxFieldName = getUniqueFieldName( mAttributePrefix +
"max" );
164 QgsField maxField( maxFieldName, QVariant::Double,
"double precision" );
170 rangeFieldName = getUniqueFieldName( mAttributePrefix +
"range" );
171 QgsField rangeField( rangeFieldName, QVariant::Double,
"double precision" );
177 minorityFieldName = getUniqueFieldName( mAttributePrefix +
"minority" );
178 QgsField minorityField( minorityFieldName, QVariant::Double,
"double precision" );
184 majorityFieldName = getUniqueFieldName( mAttributePrefix +
"majority" );
185 QgsField majField( majorityFieldName, QVariant::Double,
"double precision" );
191 varietyFieldName = getUniqueFieldName( mAttributePrefix +
"variety" );
192 QgsField varietyField( varietyFieldName, QVariant::Int,
"int" );
198 int countIndex = mStatistics & QgsZonalStatistics::Count ? vectorProvider->
fieldNameIndex( countFieldName ) : -1;
199 int sumIndex = mStatistics & QgsZonalStatistics::Sum ? vectorProvider->
fieldNameIndex( sumFieldName ) : -1;
200 int meanIndex = mStatistics & QgsZonalStatistics::Mean ? vectorProvider->
fieldNameIndex( meanFieldName ) : -1;
201 int medianIndex = mStatistics & QgsZonalStatistics::Median ? vectorProvider->
fieldNameIndex( medianFieldName ) : -1;
202 int stdevIndex = mStatistics & QgsZonalStatistics::StDev ? vectorProvider->
fieldNameIndex( stdevFieldName ) : -1;
203 int minIndex = mStatistics & QgsZonalStatistics::Min ? vectorProvider->
fieldNameIndex( minFieldName ) : -1;
204 int maxIndex = mStatistics & QgsZonalStatistics::Max ? vectorProvider->
fieldNameIndex( maxFieldName ) : -1;
205 int rangeIndex = mStatistics & QgsZonalStatistics::Range ? vectorProvider->
fieldNameIndex( rangeFieldName ) : -1;
206 int minorityIndex = mStatistics & QgsZonalStatistics::Minority ? vectorProvider->
fieldNameIndex( minorityFieldName ) : -1;
207 int majorityIndex = mStatistics & QgsZonalStatistics::Majority ? vectorProvider->
fieldNameIndex( majorityFieldName ) : -1;
208 int varietyIndex = mStatistics & QgsZonalStatistics::Variety ? vectorProvider->
fieldNameIndex( varietyFieldName ) : -1;
210 if (( mStatistics & QgsZonalStatistics::Count && countIndex == -1 )
211 || ( mStatistics & QgsZonalStatistics::Sum && sumIndex == -1 )
212 || ( mStatistics & QgsZonalStatistics::Mean && meanIndex == -1 )
213 || ( mStatistics & QgsZonalStatistics::Median && medianIndex == -1 )
214 || ( mStatistics & QgsZonalStatistics::StDev && stdevIndex == -1 )
215 || ( mStatistics & QgsZonalStatistics::Min && minIndex == -1 )
216 || ( mStatistics & QgsZonalStatistics::Max && maxIndex == -1 )
217 || ( mStatistics & QgsZonalStatistics::Range && rangeIndex == -1 )
218 || ( mStatistics & QgsZonalStatistics::Minority && minorityIndex == -1 )
219 || ( mStatistics & QgsZonalStatistics::Majority && majorityIndex == -1 )
220 || ( mStatistics & QgsZonalStatistics::Variety && varietyIndex == -1 )
242 ( mStatistics & QgsZonalStatistics::StDev );
244 ( mStatistics & QgsZonalStatistics::Majority );
246 FeatureStats featureStats( statsStoreValues, statsStoreValueCount );
247 int featureCounter = 0;
276 int offsetX, offsetY, nCellsX, nCellsY;
277 if ( cellInfoForBBox( rasterBBox, featureRect, cellsizeX, cellsizeY, offsetX, offsetY, nCellsX, nCellsY ) != 0 )
284 if (( offsetX + nCellsX ) > nCellsXGDAL )
286 nCellsX = nCellsXGDAL - offsetX;
288 if (( offsetY + nCellsY ) > nCellsYGDAL )
290 nCellsY = nCellsYGDAL - offsetY;
293 statisticsFromMiddlePointTest( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
294 rasterBBox, featureStats );
296 if ( featureStats.count <= 1 )
299 statisticsFromPreciseIntersection( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
300 rasterBBox, featureStats );
305 if ( mStatistics & QgsZonalStatistics::Count )
306 changeAttributeMap.
insert( countIndex,
QVariant( featureStats.count ) );
307 if ( mStatistics & QgsZonalStatistics::Sum )
308 changeAttributeMap.
insert( sumIndex,
QVariant( featureStats.sum ) );
309 if ( featureStats.count > 0 )
311 double mean = featureStats.sum / featureStats.count;
312 if ( mStatistics & QgsZonalStatistics::Mean )
314 if ( mStatistics & QgsZonalStatistics::Median )
316 qSort( featureStats.values.begin(), featureStats.values.end() );
317 int size = featureStats.values.count();
318 bool even = ( size % 2 ) < 1;
322 medianValue = ( featureStats.values[size / 2 - 1] + featureStats.values[size / 2] ) / 2;
326 medianValue = featureStats.values[( size + 1 ) / 2 - 1];
328 changeAttributeMap.
insert( medianIndex,
QVariant( medianValue ) );
330 if ( mStatistics & QgsZonalStatistics::StDev )
332 double sumSquared = 0;
333 for (
int i = 0; i < featureStats.values.count(); ++i )
335 double diff = featureStats.values.at( i ) - mean;
336 sumSquared += diff * diff;
338 double stdev = qPow( sumSquared / featureStats.values.count(), 0.5 );
341 if ( mStatistics & QgsZonalStatistics::Min )
342 changeAttributeMap.
insert( minIndex,
QVariant( featureStats.min ) );
343 if ( mStatistics & QgsZonalStatistics::Max )
344 changeAttributeMap.
insert( maxIndex,
QVariant( featureStats.max ) );
345 if ( mStatistics & QgsZonalStatistics::Range )
346 changeAttributeMap.
insert( rangeIndex,
QVariant( featureStats.max - featureStats.min ) );
347 if ( mStatistics & QgsZonalStatistics::Minority || mStatistics & QgsZonalStatistics::Majority )
349 QList<int> vals = featureStats.valueCount.values();
351 if ( mStatistics & QgsZonalStatistics::Minority )
353 float minorityKey = featureStats.valueCount.key( vals.
first() );
354 changeAttributeMap.
insert( minorityIndex,
QVariant( minorityKey ) );
356 if ( mStatistics & QgsZonalStatistics::Majority )
358 float majKey = featureStats.valueCount.key( vals.
last() );
362 if ( mStatistics & QgsZonalStatistics::Variety )
363 changeAttributeMap.
insert( varietyIndex,
QVariant( featureStats.valueCount.count() ) );
366 changeMap.
insert( f.
id(), changeAttributeMap );
377 GDALClose( inputDataset );
388 int QgsZonalStatistics::cellInfoForBBox(
const QgsRectangle& rasterBBox,
const QgsRectangle& featureBBox,
double cellSizeX,
double cellSizeY,
389 int& offsetX,
int& offsetY,
int& nCellsX,
int& nCellsY )
const
395 nCellsX = 0; nCellsY = 0; offsetX = 0; offsetY = 0;
400 offsetX = ( int )(( intersectBox.
xMinimum() - rasterBBox.
xMinimum() ) / cellSizeX );
401 offsetY = ( int )(( rasterBBox.
yMaximum() - intersectBox.
yMaximum() ) / cellSizeY );
403 int maxColumn = ( int )(( intersectBox.
xMaximum() - rasterBBox.
xMinimum() ) / cellSizeX ) + 1;
404 int maxRow = ( int )(( rasterBBox.
yMaximum() - intersectBox.
yMinimum() ) / cellSizeY ) + 1;
406 nCellsX = maxColumn - offsetX;
407 nCellsY = maxRow - offsetY;
412 void QgsZonalStatistics::statisticsFromMiddlePointTest(
void* band,
const QgsGeometry* poly,
int pixelOffsetX,
413 int pixelOffsetY,
int nCellsX,
int nCellsY,
double cellSizeX,
double cellSizeY,
const QgsRectangle& rasterBBox, FeatureStats &stats )
415 double cellCenterX, cellCenterY;
417 float* scanLine = (
float * ) CPLMalloc(
sizeof(
float ) * nCellsX );
418 cellCenterY = rasterBBox.
yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
421 const GEOSGeometry* polyGeos = poly->
asGeos();
428 const GEOSPreparedGeometry* polyGeosPrepared = GEOSPrepare_r( geosctxt, poly->
asGeos() );
429 if ( !polyGeosPrepared )
434 GEOSCoordSequence* cellCenterCoords = 0;
435 GEOSGeometry* currentCellCenter = 0;
437 for (
int i = 0; i < nCellsY; ++i )
439 if ( GDALRasterIO( band, GF_Read, pixelOffsetX, pixelOffsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 )
444 cellCenterX = rasterBBox.
xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;
445 for (
int j = 0; j < nCellsX; ++j )
447 if ( validPixel( scanLine[j] ) )
449 GEOSGeom_destroy_r( geosctxt, currentCellCenter );
450 cellCenterCoords = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
451 GEOSCoordSeq_setX_r( geosctxt, cellCenterCoords, 0, cellCenterX );
452 GEOSCoordSeq_setY_r( geosctxt, cellCenterCoords, 0, cellCenterY );
453 currentCellCenter = GEOSGeom_createPoint_r( geosctxt, cellCenterCoords );
454 if ( GEOSPreparedContains_r( geosctxt, polyGeosPrepared, currentCellCenter ) )
456 stats.addValue( scanLine[j] );
459 cellCenterX += cellSizeX;
461 cellCenterY -= cellSizeY;
463 GEOSGeom_destroy_r( geosctxt, currentCellCenter );
465 GEOSPreparedGeom_destroy_r( geosctxt, polyGeosPrepared );
468 void QgsZonalStatistics::statisticsFromPreciseIntersection(
void* band,
const QgsGeometry* poly,
int pixelOffsetX,
469 int pixelOffsetY,
int nCellsX,
int nCellsY,
double cellSizeX,
double cellSizeY,
const QgsRectangle& rasterBBox, FeatureStats &stats )
473 double currentY = rasterBBox.
yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
474 float* pixelData = (
float * ) CPLMalloc(
sizeof(
float ) );
477 double hCellSizeX = cellSizeX / 2.0;
478 double hCellSizeY = cellSizeY / 2.0;
479 double pixelArea = cellSizeX * cellSizeY;
482 for (
int row = 0; row < nCellsY; ++row )
484 double currentX = rasterBBox.
xMinimum() + cellSizeX / 2.0 + pixelOffsetX * cellSizeX;
485 for (
int col = 0; col < nCellsX; ++col )
487 GDALRasterIO( band, GF_Read, pixelOffsetX + col, pixelOffsetY + row, nCellsX, 1, pixelData, 1, 1, GDT_Float32, 0, 0 );
488 if ( !validPixel( *pixelData ) )
492 if ( pixelRectGeometry )
496 if ( intersectGeometry )
498 double intersectionArea = intersectGeometry->
area();
499 if ( intersectionArea >= 0.0 )
501 weight = intersectionArea / pixelArea;
502 stats.addValue( *pixelData, weight );
504 delete intersectGeometry;
506 delete pixelRectGeometry;
507 pixelRectGeometry = 0;
509 currentX += cellSizeX;
511 currentY -= cellSizeY;
513 CPLFree( pixelData );
516 bool QgsZonalStatistics::validPixel(
float value )
const
518 if ( value == mInputNodataValue || qIsNaN( value ) )
525 QString QgsZonalStatistics::getUniqueFieldName(
QString fieldName )
538 for (
int idx = 0; idx < providerFields.
count(); ++idx )
540 if ( shortName == providerFields[idx].name() )
553 shortName =
QString(
"%1_%2" ).
arg( fieldName.
mid( 0, 8 ) ).arg( n );
558 for (
int idx = 0; idx < providerFields.
count(); ++idx )
560 if ( shortName == providerFields[idx].name() )
565 shortName =
QString(
"%1_%2" ).
arg( fieldName.
mid( 0, 8 ) ).arg( n );
569 shortName =
QString(
"%1_%2" ).
arg( fieldName.
mid( 0, 7 ) ).arg( n );
QgsFeatureId id() const
Get the feature ID for this feature.
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.
bool isEmpty() const
test if rectangle is empty.
virtual bool addAttributes(const QList< QgsField > &attributes)
Adds new attributes.
void setMaximum(int maximum)
void push_back(const T &value)
double yMaximum() const
Get the y maximum value (top side of rectangle)
QgsRectangle boundingBox() const
Returns the bounding box of this feature.
QgsFeatureRequest & setSubsetOfAttributes(const QgsAttributeList &attrs)
Set a subset of attributes that will be fetched.
Container of fields for a vector layer.
A geometry is the spatial representation of a feature.
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
int fieldNameIndex(const QString &fieldName) const
Returns the index of a field name or -1 if the field does not exist.
const GEOSGeometry * asGeos() const
Returns a geos geometry.
Minority of pixel values.
Variety (count of distinct) pixel values.
void setValue(int progress)
double yMinimum() const
Get the y minimum value (bottom side of rectangle)
double xMaximum() const
Get the x maximum value (right side of rectangle)
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.
Majority of pixel values.
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.
QGis::GeometryType geometryType() const
Returns point, line or polygon.
Encapsulate a field in an attribute table or data source.
QgsZonalStatistics(QgsVectorLayer *polygonLayer, const QString &rasterFile, const QString &attributePrefix="", int rasterBand=1, Statistics stats=Statistics(Count|Sum|Mean))
bool contains(QChar ch, Qt::CaseSensitivity cs) const
QgsGeometry * intersection(const QgsGeometry *geometry) const
Returns a geometry representing the points shared by this geometry and other.
Range of pixel values (max - min)
virtual const QgsFields & fields() const =0
Return a map of indexes with field names for this layer.
int calculateStatistics(QProgressDialog *p)
Starts the calculation.
QString mid(int position, int n) const
static GEOSContextHandle_t getGEOSHandler()
return GEOS context handle
QgsRectangle intersect(const QgsRectangle *rect) const
return the intersection with the given rectangle
const QgsGeometry * constGeometry() const
Gets a const pointer to the geometry object associated with this feature.
static QgsGeometry * fromRect(const QgsRectangle &rect)
Creates a new geometry from a QgsRectangle.
A class that calculates raster statistics (count, sum, mean) for a polygon or multipolygon layer and ...
iterator insert(const Key &key, const T &value)
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.
QString arg(qlonglong a, int fieldWidth, int base, const QChar &fillChar) const
double xMinimum() const
Get the x minimum value (left side of rectangle)
double area() const
Returns the area of the geometry using GEOS.
Standard deviation of pixel values.