[QGIS Commit] r11931 - in trunk/qgis/src/analysis: . vector
svn_qgis at osgeo.org
svn_qgis at osgeo.org
Fri Nov 6 06:23:18 EST 2009
Author: mhugent
Date: 2009-11-06 06:23:17 -0500 (Fri, 06 Nov 2009)
New Revision: 11931
Added:
trunk/qgis/src/analysis/vector/qgszonalstatistics.cpp
trunk/qgis/src/analysis/vector/qgszonalstatistics.h
Modified:
trunk/qgis/src/analysis/CMakeLists.txt
Log:
Added zonal statistics class to analysis lib
Modified: trunk/qgis/src/analysis/CMakeLists.txt
===================================================================
--- trunk/qgis/src/analysis/CMakeLists.txt 2009-11-06 10:55:35 UTC (rev 11930)
+++ trunk/qgis/src/analysis/CMakeLists.txt 2009-11-06 11:23:17 UTC (rev 11931)
@@ -25,6 +25,7 @@
raster/qgsaspectfilter.cpp
raster/qgstotalcurvaturefilter.cpp
vector/qgsgeometryanalyzer.cpp
+ vector/qgszonalstatistics.cpp
)
SET(QGIS_ANALYSIS_MOC_HDRS
@@ -93,7 +94,7 @@
# Added by Tim to install headers
-SET(QGIS_ANALYSIS_HDRS vector/qgsgeometryanalyzer.h
+SET(QGIS_ANALYSIS_HDRS vector/qgsgeometryanalyzer.h vector/qgszonalstatistics.h
)
INSTALL(CODE "MESSAGE(\"Installing ANALYSIS headers...\")")
Added: trunk/qgis/src/analysis/vector/qgszonalstatistics.cpp
===================================================================
--- trunk/qgis/src/analysis/vector/qgszonalstatistics.cpp (rev 0)
+++ trunk/qgis/src/analysis/vector/qgszonalstatistics.cpp 2009-11-06 11:23:17 UTC (rev 11931)
@@ -0,0 +1,252 @@
+/***************************************************************************
+ qgszonalstatistics.cpp - description
+ ----------------------------
+ begin : August 29th, 2009
+ copyright : (C) 2009 by Marco Hugentobler
+ email : marco at hugis dot net
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ ***************************************************************************/
+
+#include "qgszonalstatistics.h"
+#include "qgsgeometry.h"
+#include "qgsvectordataprovider.h"
+#include "qgsvectorlayer.h"
+#include "gdal.h"
+#include "cpl_string.h"
+#include <QProgressDialog>
+
+QgsZonalStatistics::QgsZonalStatistics(QgsVectorLayer* polygonLayer, const QString& rasterFile, const QString& attributePrefix, int rasterBand): \
+ mRasterFilePath(rasterFile), mRasterBand(rasterBand), mPolygonLayer(polygonLayer), mAttributePrefix(attributePrefix), mInputNodataValue( -1 )
+{
+
+}
+
+QgsZonalStatistics::QgsZonalStatistics(): mRasterBand(0), mPolygonLayer(0)
+{
+
+}
+
+QgsZonalStatistics::~QgsZonalStatistics()
+{
+
+}
+
+int QgsZonalStatistics::calculateStatistics(QProgressDialog* p)
+{
+ if(!mPolygonLayer || !mPolygonLayer->geometryType() == QGis::Polygon)
+ {
+ return 1;
+ }
+
+ QgsVectorDataProvider* vectorProvider = mPolygonLayer->dataProvider();
+ if(!vectorProvider)
+ {
+ return 2;
+ }
+
+ //open the raster layer and the raster band
+ GDALAllRegister();
+ GDALDatasetH inputDataset = GDALOpen(mRasterFilePath.toLocal8Bit().data(), GA_ReadOnly );
+ if ( inputDataset == NULL )
+ {
+ return 3;
+ }
+
+ if ( GDALGetRasterCount( inputDataset ) < (mRasterBand - 1) )
+ {
+ GDALClose( inputDataset );
+ return 4;
+ }
+
+ GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset, mRasterBand );
+ if ( rasterBand == NULL )
+ {
+ GDALClose( inputDataset );
+ return 5;
+ }
+ mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, NULL );
+
+ //get geometry info about raster layer
+ int nCellsX = GDALGetRasterXSize( inputDataset );
+ int nCellsY = GDALGetRasterYSize( inputDataset );
+ double geoTransform[6];
+ if ( GDALGetGeoTransform( inputDataset, geoTransform ) != CE_None )
+ {
+ GDALClose( inputDataset );
+ return 6;
+ }
+ double cellsizeX = geoTransform[1];
+ if(cellsizeX < 0)
+ {
+ cellsizeX = -cellsizeX;
+ }
+ double cellsizeY = geoTransform[5];
+ if(cellsizeY < 0)
+ {
+ cellsizeY = -cellsizeY;
+ }
+ QgsRectangle rasterBBox(geoTransform[0], geoTransform[3] - (nCellsY * cellsizeY), geoTransform[0] + (nCellsX * cellsizeX), geoTransform[3]);
+
+ //add the new count, sum, mean fields to the provider
+ QList<QgsField> newFieldList;
+ QgsField countField(mAttributePrefix + "count", QVariant::Int);
+ QgsField sumField(mAttributePrefix + "sum", QVariant::Double);
+ QgsField meanField(mAttributePrefix + "mean", QVariant::Double);
+ newFieldList.push_back(countField);
+ newFieldList.push_back(sumField);
+ newFieldList.push_back(meanField);
+ if(!vectorProvider->addAttributes( newFieldList ))
+ {
+ return 7;
+ }
+
+ //index of the new fields
+ int countIndex = vectorProvider->fieldNameIndex(mAttributePrefix + "count");
+ int sumIndex = vectorProvider->fieldNameIndex(mAttributePrefix + "sum");
+ int meanIndex = vectorProvider->fieldNameIndex(mAttributePrefix + "mean");
+
+ if (countIndex == -1 || sumIndex == -1 || meanIndex == -1)
+ {
+ return 8;
+ }
+
+ //progress dialog
+ long featureCount = vectorProvider->featureCount();
+ if(p)
+ {
+ p->setMaximum(featureCount);
+ }
+
+
+ //iterate over each polygon
+ vectorProvider->select(QgsAttributeList(), QgsRectangle(), true, false);
+ vectorProvider->rewind();
+ QgsFeature f;
+ double count = 0;
+ double sum = 0;
+ double mean = 0;
+ float* scanLine;
+ int featureCounter = 0;
+ //x- and y- coordinate of current cell
+ double cellCenterX = 0;
+ double cellCenterY = 0;
+ QgsPoint currentCellCenter;
+
+ while(vectorProvider->nextFeature(f))
+ {
+ if(p)
+ {
+ p->setValue(featureCounter);
+ }
+
+ if(p && p->wasCanceled())
+ {
+ break;
+ }
+
+ QgsGeometry* featureGeometry = f.geometry();
+ if(!featureGeometry)
+ {
+ ++featureCounter;
+ continue;
+ }
+
+ int offsetX, offsetY, nCellsX, nCellsY;
+ if(cellInfoForBBox(rasterBBox, featureGeometry->boundingBox(), cellsizeX, cellsizeY, offsetX, offsetY, nCellsX, nCellsY) != 0)
+ {
+ ++featureCounter;
+ continue;
+ }
+
+ scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );
+ cellCenterY = rasterBBox.yMaximum() - offsetY * cellsizeY - cellsizeY / 2;
+ count = 0;
+ sum = 0;
+ float currentValue;
+
+ for(int i = 0; i < nCellsY; ++i)
+ {
+ GDALRasterIO( rasterBand, GF_Read, offsetX, offsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 );
+ cellCenterX = rasterBBox.xMinimum() + offsetX * cellsizeX + cellsizeX / 2;
+ for(int j = 0; j < nCellsX; ++j)
+ {
+ currentCellCenter = QgsPoint(cellCenterX, cellCenterY);
+ if(featureGeometry->contains(¤tCellCenter))
+ {
+ if(scanLine[j] != mInputNodataValue) //don't consider nodata values
+ {
+ sum += scanLine[j];
+ ++count;
+ }
+ }
+ cellCenterX += cellsizeX;
+ }
+ cellCenterY -= cellsizeY;
+ }
+
+ if(count == 0)
+ {
+ mean = 0;
+ }
+ else
+ {
+ mean = sum / count;
+ }
+
+ //write the new AEY value to the vector data provider
+ QgsChangedAttributesMap changeMap;
+ QgsAttributeMap changeAttributeMap;
+ changeAttributeMap.insert(countIndex, QVariant(count));
+ changeAttributeMap.insert(sumIndex, QVariant(sum));
+ changeAttributeMap.insert(meanIndex, QVariant(mean));
+ changeMap.insert(f.id(), changeAttributeMap);
+ vectorProvider->changeAttributeValues(changeMap);
+
+ CPLFree(scanLine);
+ ++featureCounter;
+ }
+
+ if(p)
+ {
+ p->setValue(featureCount);
+ }
+
+ GDALClose( inputDataset );
+ return 0;
+}
+
+int QgsZonalStatistics::cellInfoForBBox(const QgsRectangle& rasterBBox, const QgsRectangle& featureBBox, double cellSizeX, double cellSizeY, \
+ int& offsetX, int& offsetY, int& nCellsX, int& nCellsY) const
+{
+ //get intersecting bbox
+ QgsRectangle intersectBox = rasterBBox.intersect(&featureBBox);
+ if(intersectBox.isEmpty())
+ {
+ nCellsX = 0; nCellsY = 0; offsetX = 0; offsetY = 0;
+ return 0;
+ }
+
+ //get offset in pixels in x- and y- direction
+ offsetX = (int)( (intersectBox.xMinimum() - rasterBBox.xMinimum()) / cellSizeX);
+ offsetY = (int)( (rasterBBox.yMaximum() - intersectBox.yMaximum()) / cellSizeY);
+
+ int maxColumn = (int) ( (intersectBox.xMaximum() - rasterBBox.xMinimum()) / cellSizeX) + 1;
+ int maxRow = (int) ( (rasterBBox.yMaximum() - intersectBox.yMinimum()) / cellSizeY ) + 1;
+
+ nCellsX = maxColumn - offsetX;
+ nCellsY = maxRow - offsetY;
+
+ return 0;
+}
+
+
+
+
Added: trunk/qgis/src/analysis/vector/qgszonalstatistics.h
===================================================================
--- trunk/qgis/src/analysis/vector/qgszonalstatistics.h (rev 0)
+++ trunk/qgis/src/analysis/vector/qgszonalstatistics.h 2009-11-06 11:23:17 UTC (rev 11931)
@@ -0,0 +1,54 @@
+/***************************************************************************
+ qgszonalstatistics.h - description
+ ----------------------------
+ begin : August 29th, 2009
+ copyright : (C) 2009 by Marco Hugentobler
+ email : marco at hugis dot net
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ ***************************************************************************/
+
+#ifndef QGSZONALSTATISTICS_H
+#define QGSZONALSTATISTICS_H
+
+#include "qgsrectangle.h"
+#include <QString>
+
+class QgsVectorLayer;
+class QProgressDialog;
+
+/**A class that calculates raster statistics (count, sum, mean) for a polygon or multipolygon layer and appends the results as attributes*/
+class ANALYSIS_EXPORT QgsZonalStatistics
+{
+ public:
+ QgsZonalStatistics(QgsVectorLayer* polygonLayer, const QString& rasterFile, const QString& attributePrefix = "", int rasterBand = 1);
+ ~QgsZonalStatistics();
+
+ /**Starts the calculation
+ @return 0 in case of success*/
+ int calculateStatistics(QProgressDialog* p);
+
+ private:
+ QgsZonalStatistics();
+ /**Analysis what cells need to be considered to cover the bounding box of a feature
+ @return 0 in case of success*/
+ int cellInfoForBBox(const QgsRectangle& rasterBBox, const QgsRectangle& featureBBox, double cellSizeX, double cellSizeY, \
+ int& offsetX, int& offsetY, int& nCellsX, int& nCellsY) const;
+
+ QString mRasterFilePath;
+ /**Raster band to calculate statistics from (defaults to 1)*/
+ int mRasterBand;
+ QgsVectorLayer* mPolygonLayer;
+ QString mAttributePrefix;
+ /**The nodata value of the input layer*/
+ float mInputNodataValue;
+};
+
+#endif // QGSZONALSTATISTICS_H
More information about the QGIS-commit
mailing list