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.
qgis · Jun 8, 2018 · 985aee6 · 985aee6
1 parent e9bf7f1
commit 985aee6
Show file tree

Hide file tree

Showing 3 changed files with 33 additions and 36 deletions.
diff --git a/src/analysis/vector/qgszonalstatistics.cpp b/src/analysis/vector/qgszonalstatistics.cpp
@@ -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 );

diff --git a/src/analysis/vector/qgszonalstatistics.h b/src/analysis/vector/qgszonalstatistics.h
@@ -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

diff --git a/tests/src/analysis/testqgszonalstatistics.cpp b/tests/src/analysis/testqgszonalstatistics.cpp
@@ -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 );