Skip to content
Permalink
Browse files

Fix incorrect calculation of zonal statistics

Raster block extent was based on intersection of raster vs
feature's bounding box, which was not necessarily snapped to
multiples of the pixel size. But the pixel center point/extent
was being calculated as though the retrieved extent was an
exact multiple of the pixel size. This led to incorrect
retrieval of pixel values.
  • Loading branch information
nyalldawson committed Jun 8, 2018
1 parent e9bf7f1 commit 985aee6920bc0371bc299e576c750d8b6b15d300
@@ -258,17 +258,18 @@ int QgsZonalStatistics::calculateStatistics( QgsFeedback *feedback )
continue;
}

int offsetX, offsetY, nCellsX, nCellsY;
cellInfoForBBox( rasterBBox, featureRect, mCellSizeX, mCellSizeY, offsetX, offsetY, nCellsX, nCellsY, nCellsXProvider, nCellsYProvider );
int nCellsX, nCellsY;
QgsRectangle rasterBlockExtent;
cellInfoForBBox( rasterBBox, featureRect, mCellSizeX, mCellSizeY, nCellsX, nCellsY, nCellsXProvider, nCellsYProvider, rasterBlockExtent );

statisticsFromMiddlePointTest( featureGeometry, offsetX, offsetY, nCellsX, nCellsY, mCellSizeX, mCellSizeY,
rasterBBox, featureStats );
statisticsFromMiddlePointTest( featureGeometry, nCellsX, nCellsY, mCellSizeX, mCellSizeY,
rasterBlockExtent, featureStats );

if ( featureStats.count <= 1 )
{
//the cell resolution is probably larger than the polygon area. We switch to precise pixel - polygon intersection in this case
statisticsFromPreciseIntersection( featureGeometry, offsetX, offsetY, nCellsX, nCellsY, mCellSizeX, mCellSizeY,
rasterBBox, featureStats );
statisticsFromPreciseIntersection( featureGeometry, nCellsX, nCellsY, mCellSizeX, mCellSizeY,
rasterBlockExtent, featureStats );
}

//write the statistics value to the vector data provider
@@ -362,22 +363,21 @@ int QgsZonalStatistics::calculateStatistics( QgsFeedback *feedback )
}

void QgsZonalStatistics::cellInfoForBBox( const QgsRectangle &rasterBBox, const QgsRectangle &featureBBox, double cellSizeX, double cellSizeY,
int &offsetX, int &offsetY, int &nCellsX, int &nCellsY, int rasterWidth, int rasterHeight ) const
int &nCellsX, int &nCellsY, int rasterWidth, int rasterHeight, QgsRectangle &rasterBlockExtent ) const
{
//get intersecting bbox
QgsRectangle intersectBox = rasterBBox.intersect( &featureBBox );
if ( intersectBox.isEmpty() )
{
nCellsX = 0;
nCellsY = 0;
offsetX = 0;
offsetY = 0;
rasterBlockExtent = QgsRectangle();
return;
}

//get offset in pixels in x- and y- direction
offsetX = static_cast< int >( std::floor( ( intersectBox.xMinimum() - rasterBBox.xMinimum() ) / cellSizeX ) );
offsetY = static_cast< int >( std::floor( ( rasterBBox.yMaximum() - intersectBox.yMaximum() ) / cellSizeY ) );
int offsetX = static_cast< int >( std::floor( ( intersectBox.xMinimum() - rasterBBox.xMinimum() ) / cellSizeX ) );
int offsetY = static_cast< int >( std::floor( ( rasterBBox.yMaximum() - intersectBox.yMaximum() ) / cellSizeY ) );

int maxColumn = static_cast< int >( std::floor( ( intersectBox.xMaximum() - rasterBBox.xMinimum() ) / cellSizeX ) ) + 1;
int maxRow = static_cast< int >( std::floor( ( rasterBBox.yMaximum() - intersectBox.yMinimum() ) / cellSizeY ) ) + 1;
@@ -388,14 +388,18 @@ void QgsZonalStatistics::cellInfoForBBox( const QgsRectangle &rasterBBox, const
//avoid access to cells outside of the raster (may occur because of rounding)
nCellsX = std::min( offsetX + nCellsX, rasterWidth ) - offsetX;
nCellsY = std::min( offsetY + nCellsY, rasterHeight ) - offsetY;

rasterBlockExtent = QgsRectangle( rasterBBox.xMinimum() + offsetX * cellSizeX,
rasterBBox.yMaximum() - offsetY * cellSizeY,
rasterBBox.xMinimum() + ( nCellsX + offsetX ) * cellSizeX,
rasterBBox.yMaximum() - ( nCellsY + offsetY ) * cellSizeY );
}

void QgsZonalStatistics::statisticsFromMiddlePointTest( const QgsGeometry &poly, int pixelOffsetX,
int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle &rasterBBox, FeatureStats &stats )
void QgsZonalStatistics::statisticsFromMiddlePointTest( const QgsGeometry &poly, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle &rasterBBox, FeatureStats &stats )
{
double cellCenterX, cellCenterY;

cellCenterY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
cellCenterY = rasterBBox.yMaximum() - cellSizeY / 2;
stats.reset();

std::unique_ptr< QgsGeometryEngine > polyEngine( QgsGeometry::createGeometryEngine( poly.constGet( ) ) );
@@ -405,13 +409,10 @@ void QgsZonalStatistics::statisticsFromMiddlePointTest( const QgsGeometry &poly,
}
polyEngine->prepareGeometry();

QgsRectangle featureBBox = poly.boundingBox().intersect( &rasterBBox );
QgsRectangle intersectBBox = rasterBBox.intersect( &featureBBox );

std::unique_ptr< QgsRasterBlock > block( mRasterInterface->block( mRasterBand, intersectBBox, nCellsX, nCellsY ) );
std::unique_ptr< QgsRasterBlock > block( mRasterInterface->block( mRasterBand, rasterBBox, nCellsX, nCellsY ) );
for ( int i = 0; i < nCellsY; ++i )
{
cellCenterX = rasterBBox.xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;
cellCenterX = rasterBBox.xMinimum() + cellSizeX / 2;
for ( int j = 0; j < nCellsX; ++j )
{
double pixelValue = block->value( i, j );
@@ -429,33 +430,29 @@ void QgsZonalStatistics::statisticsFromMiddlePointTest( const QgsGeometry &poly,
}
}

void QgsZonalStatistics::statisticsFromPreciseIntersection( const QgsGeometry &poly, int pixelOffsetX,
int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle &rasterBBox, FeatureStats &stats )
void QgsZonalStatistics::statisticsFromPreciseIntersection( const QgsGeometry &poly, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle &rasterBBox, FeatureStats &stats )
{
stats.reset();

double currentY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
double currentY = rasterBBox.yMaximum() - cellSizeY / 2;
QgsGeometry pixelRectGeometry;

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

QgsRectangle featureBBox = poly.boundingBox().intersect( &rasterBBox );
QgsRectangle intersectBBox = rasterBBox.intersect( &featureBBox );

std::unique_ptr< QgsGeometryEngine > polyEngine( QgsGeometry::createGeometryEngine( poly.constGet( ) ) );
if ( !polyEngine )
{
return;
}
polyEngine->prepareGeometry();

std::unique_ptr< QgsRasterBlock > block( mRasterInterface->block( mRasterBand, intersectBBox, nCellsX, nCellsY ) );
std::unique_ptr< QgsRasterBlock > block( mRasterInterface->block( mRasterBand, rasterBBox, nCellsX, nCellsY ) );
for ( int i = 0; i < nCellsY; ++i )
{
double currentX = rasterBBox.xMinimum() + cellSizeX / 2.0 + pixelOffsetX * cellSizeX;
double currentX = rasterBBox.xMinimum() + cellSizeX / 2.0;
for ( int j = 0; j < nCellsX; ++j )
{
double pixelValue = block->value( i, j );
@@ -155,16 +155,16 @@ class ANALYSIS_EXPORT QgsZonalStatistics
/**
* Analyzes which cells need to be considered to completely cover the bounding box of a feature.
*/
void cellInfoForBBox( const QgsRectangle &rasterBBox, const QgsRectangle &featureBBox, double cellSizeX, double cellSizeY,
int &offsetX, int &offsetY, int &nCellsX, int &nCellsY,
int rasterWidth, int rasterHeight ) const;
void cellInfoForBBox( const QgsRectangle &rasterBBox, const QgsRectangle &featureBBox, double cellSizeX, double cellSizeY, int &nCellsX, int &nCellsY,
int rasterWidth, int rasterHeight,
QgsRectangle &rasterBlockExtent ) const;

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

//! Returns statistics with precise pixel - polygon intersection test (slow)
void statisticsFromPreciseIntersection( const QgsGeometry &poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY,
void statisticsFromPreciseIntersection( const QgsGeometry &poly, int nCellsX, int nCellsY,
double cellSizeX, double cellSizeY, const QgsRectangle &rasterBBox, FeatureStats &stats );

//! Tests whether a pixel's value should be included in the result
@@ -270,8 +270,8 @@ void TestQgsZonalStatistics::testNoData()

fetched = it.nextFeature( f );
QVERIFY( fetched );
QCOMPARE( f.attribute( "ncount" ).toDouble(), 103.0 );
QCOMPARE( f.attribute( "nsum" ).toDouble(), 90536.0 );
QCOMPARE( f.attribute( "ncount" ).toDouble(), 50.0 );
QCOMPARE( f.attribute( "nsum" ).toDouble(), 43868.0 );

fetched = it.nextFeature( f );
QVERIFY( fetched );
@@ -293,8 +293,8 @@ void TestQgsZonalStatistics::testNoData()

fetched = it.nextFeature( f );
QVERIFY( fetched );
QCOMPARE( f.attribute( "uncount" ).toDouble(), 52.0 );
QCOMPARE( f.attribute( "unsum" ).toDouble(), 45374.0 );
QCOMPARE( f.attribute( "uncount" ).toDouble(), 25.0 );
QCOMPARE( f.attribute( "unsum" ).toDouble(), 21750.0 );

fetched = it.nextFeature( f );
QVERIFY( fetched );

0 comments on commit 985aee6

Please sign in to comment.
You can’t perform that action at this time.