Fix calculation of zonal stats when source contains nodata or nan

pixels (fix #11135)

Author: nyalldawson

src/analysis/vector/qgszonalstatistics.cpp

@@ -302,27 +302,24 @@ void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, QgsGeometry*
     cellCenterX = rasterBBox.xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;
     for ( int j = 0; j < nCellsX; ++j )
     {
-      GEOSGeom_destroy_r( geosctxt, currentCellCenter );
-      cellCenterCoords = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
-      GEOSCoordSeq_setX_r( geosctxt, cellCenterCoords, 0, cellCenterX );
-      GEOSCoordSeq_setY_r( geosctxt, cellCenterCoords, 0, cellCenterY );
-      currentCellCenter = GEOSGeom_createPoint_r( geosctxt, cellCenterCoords );
-
-      if ( scanLine[j] != mInputNodataValue ) //don't consider nodata values
+      if ( validPixel( scanLine[j] ) )
       {
+        GEOSGeom_destroy_r( geosctxt, currentCellCenter );
+        cellCenterCoords = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
+        GEOSCoordSeq_setX_r( geosctxt, cellCenterCoords, 0, cellCenterX );
+        GEOSCoordSeq_setY_r( geosctxt, cellCenterCoords, 0, cellCenterY );
+        currentCellCenter = GEOSGeom_createPoint_r( geosctxt, cellCenterCoords );
         if ( GEOSPreparedContains_r( geosctxt, polyGeosPrepared, currentCellCenter ) )
         {
-          if ( !qIsNaN( scanLine[j] ) )
-          {
-            sum += scanLine[j];
-          }
+          sum += scanLine[j];
           ++count;
         }
       }
       cellCenterX += cellSizeX;
     }
     cellCenterY -= cellSizeY;
   }
+  GEOSGeom_destroy_r( geosctxt, currentCellCenter );
   CPLFree( scanLine );
   GEOSPreparedGeom_destroy_r( geosctxt, polyGeosPrepared );
 }
@@ -347,6 +344,9 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, QgsGeome
     for ( int col = 0; col < nCellsX; ++col )
     {
       GDALRasterIO( band, GF_Read, pixelOffsetX + col, pixelOffsetY + row, nCellsX, 1, pixelData, 1, 1, GDT_Float32, 0, 0 );
+      if ( !validPixel( *pixelData ) )
+        continue;
+
       pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) );
       if ( pixelRectGeometry )
       {
@@ -373,6 +373,15 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, QgsGeome
   CPLFree( pixelData );
 }
 
+bool QgsZonalStatistics::validPixel( float value ) const
+{
+  if ( value == mInputNodataValue || qIsNaN( value ) )
+  {
+    return false;
+  }
+  return true;
+}
+
 QString QgsZonalStatistics::getUniqueFieldName( QString fieldName )
 {
   QgsVectorDataProvider* dp = mPolygonLayer->dataProvider();

src/analysis/vector/qgszonalstatistics.h

@@ -54,6 +54,9 @@ class ANALYSIS_EXPORT QgsZonalStatistics
     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 );
 
+    /**Tests whether a pixel's value should be included in the result*/
+    bool validPixel( float value ) const;
+
     QString getUniqueFieldName( QString fieldName );
 
     QString mRasterFilePath;
@@ -63,6 +66,7 @@ class ANALYSIS_EXPORT QgsZonalStatistics
     QString mAttributePrefix;
     /**The nodata value of the input layer*/
     float mInputNodataValue;
+
 };
 
 #endif // QGSZONALSTATISTICS_H