Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Improve zonal statistics function in analysis lib
git-svn-id: http://svn.osgeo.org/qgis/trunk/qgis@12953 c8812cc2-4d05-0410-92ff-de0c093fc19c
  • Loading branch information
mhugent committed Feb 17, 2010
1 parent c7517f5 commit e508f3a
Show file tree
Hide file tree
Showing 2 changed files with 97 additions and 29 deletions.
116 changes: 87 additions & 29 deletions src/analysis/vector/qgszonalstatistics.cpp
Expand Up @@ -97,7 +97,7 @@ int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )

//add the new count, sum, mean fields to the provider
QList<QgsField> newFieldList;
QgsField countField( mAttributePrefix + "count", QVariant::Int );
QgsField countField( mAttributePrefix + "count", QVariant::Double );
QgsField sumField( mAttributePrefix + "sum", QVariant::Double );
QgsField meanField( mAttributePrefix + "mean", QVariant::Double );
newFieldList.push_back( countField );
Expand Down Expand Up @@ -133,12 +133,7 @@ int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
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 ) )
{
Expand Down Expand Up @@ -166,31 +161,17 @@ int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
continue;
}

scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );
cellCenterY = rasterBBox.yMaximum() - offsetY * cellsizeY - cellsizeY / 2;
count = 0;
sum = 0;
statisticsFromMiddlePointTest( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY, \
rasterBBox, sum, count );

for ( int i = 0; i < nCellsY; ++i )
if ( count <= 1 )
{
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( &currentCellCenter ) )
{
if ( scanLine[j] != mInputNodataValue ) //don't consider nodata values
{
sum += scanLine[j];
++count;
}
}
cellCenterX += cellsizeX;
}
cellCenterY -= cellsizeY;
//the cell resolution is probably larger than the polygon area. We switch to precise pixel - polygon intersection in this case
statisticsFromPreciseIntersection( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY, \
rasterBBox, sum, count );
}


if ( count == 0 )
{
mean = 0;
Expand All @@ -200,7 +181,7 @@ int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
mean = sum / count;
}

//write the new AEY value to the vector data provider
//write the statistics value to the vector data provider
QgsChangedAttributesMap changeMap;
QgsAttributeMap changeAttributeMap;
changeAttributeMap.insert( countIndex, QVariant( count ) );
Expand All @@ -209,7 +190,6 @@ int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
changeMap.insert( f.id(), changeAttributeMap );
vectorProvider->changeAttributeValues( changeMap );

CPLFree( scanLine );
++featureCounter;
}

Expand Down Expand Up @@ -246,6 +226,84 @@ int QgsZonalStatistics::cellInfoForBBox( const QgsRectangle& rasterBBox, const Q
return 0;
}

void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, QgsGeometry* poly, int pixelOffsetX, \
int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, double& sum, double& count )
{
double cellCenterX, cellCenterY;
QgsPoint currentCellCenter;

float* scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );
cellCenterY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
count = 0;
sum = 0;

for ( int i = 0; i < nCellsY; ++i )
{
GDALRasterIO( band, GF_Read, pixelOffsetX, pixelOffsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 );
cellCenterX = rasterBBox.xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;
for ( int j = 0; j < nCellsX; ++j )
{
currentCellCenter = QgsPoint( cellCenterX, cellCenterY );
if ( poly->contains( &currentCellCenter ) )
{
if ( scanLine[j] != mInputNodataValue ) //don't consider nodata values
{
sum += scanLine[j];
++count;
}
}
cellCenterX += cellSizeX;
}
cellCenterY -= cellSizeY;
}
CPLFree( scanLine );
}

void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, QgsGeometry* poly, int pixelOffsetX, \
int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, double& sum, double& count )
{
sum = 0;
count = 0;
double currentY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
float* pixelData = ( float * ) CPLMalloc( sizeof( float ) );
QgsGeometry* pixelRectGeometry = 0;
QgsGeometry* intersectGeometry = 0;

double hCellSizeX = cellSizeX / 2.0;
double hCellSizeY = cellSizeY / 2.0;
double pixelArea = cellSizeX * cellSizeY;
double intersectionArea = 0;
double weight = 0;

for ( int row = 0; row < nCellsY; ++row )
{
double currentX = rasterBBox.xMinimum() + cellSizeX / 2.0 + pixelOffsetX * cellSizeX;
for ( int col = 0; col < nCellsX; ++col )
{
GDALRasterIO( band, GF_Read, pixelOffsetX + col, pixelOffsetY + row, nCellsX, 1, pixelData, 1, 1, GDT_Float32, 0, 0 );
pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) );
if ( pixelRectGeometry )
{
//intersecton
intersectGeometry = pixelRectGeometry->intersection( poly );
if ( intersectGeometry )
{
if ( GEOSArea( intersectGeometry->asGeos(), &intersectionArea ) )
{
weight = intersectionArea / pixelArea;
count += weight;
sum += *pixelData * weight;
}
delete intersectGeometry;
}
}
currentX += cellSizeX;
}
currentY -= cellSizeY;
}
CPLFree( pixelData );
}




10 changes: 10 additions & 0 deletions src/analysis/vector/qgszonalstatistics.h
Expand Up @@ -21,6 +21,7 @@
#include "qgsrectangle.h"
#include <QString>

class QgsGeometry;
class QgsVectorLayer;
class QProgressDialog;

Expand All @@ -42,6 +43,15 @@ class ANALYSIS_EXPORT QgsZonalStatistics
int cellInfoForBBox( const QgsRectangle& rasterBBox, const QgsRectangle& featureBBox, double cellSizeX, double cellSizeY,
int& offsetX, int& offsetY, int& nCellsX, int& nCellsY ) const;

/**Returns statistics by considering the pixels where the center point is within the polygon (fast)*/
void statisticsFromMiddlePointTest( void* band, QgsGeometry* poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY, \
double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, double& sum, double& count );

/**Returns statistics with precise pixel - polygon intersection test (slow) */
void statisticsFromPreciseIntersection( void* band, QgsGeometry* poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY, \
double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, double& sum, double& count );


QString mRasterFilePath;
/**Raster band to calculate statistics from (defaults to 1)*/
int mRasterBand;
Expand Down

0 comments on commit e508f3a

Please sign in to comment.