Ensure zonal stats respsects all user set no data values

Author: nyalldawson

src/analysis/vector/qgszonalstatistics.cpp

@@ -62,7 +62,6 @@ int QgsZonalStatistics::calculateStatistics( QgsFeedback *feedback )
   }
 
   mRasterProvider = mRasterLayer->dataProvider();
-  mInputNodataValue = mRasterProvider->sourceNoDataValue( mRasterBand );
 
   //get geometry info about raster layer
   int nCellsXProvider = mRasterProvider->xSize();
@@ -402,7 +401,7 @@ void QgsZonalStatistics::statisticsFromMiddlePointTest( const QgsGeometry &poly,
     for ( int j = 0; j < nCellsX; ++j )
     {
       double pixelValue = block->value( i, j );
-      if ( validPixel( pixelValue ) )
+      if ( validPixel( pixelValue ) && !block->isNoData( i, j ) )
       {
         QgsPoint cellCenter( cellCenterX, cellCenterY );
         if ( polyEngine->contains( &cellCenter ) )
@@ -446,7 +445,7 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( const QgsGeometry &p
     for ( int j = 0; j < nCellsX; ++j )
     {
       double pixelValue = block->value( i, j );
-      if ( !validPixel( pixelValue ) )
+      if ( !validPixel( pixelValue ) || block->isNoData( i, j ) )
       {
         continue;
       }
@@ -474,9 +473,9 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( const QgsGeometry &p
   }
 }
 
-bool QgsZonalStatistics::validPixel( float value ) const
+bool QgsZonalStatistics::validPixel( double value ) const
 {
-  return !( value == mInputNodataValue || std::isnan( value ) );
+  return !std::isnan( value );
 }
 
 QString QgsZonalStatistics::getUniqueFieldName( const QString &fieldName, const QList<QgsField> &newFields )

src/analysis/vector/qgszonalstatistics.h

@@ -143,7 +143,7 @@ class ANALYSIS_EXPORT QgsZonalStatistics
                                             double cellSizeX, double cellSizeY, const QgsRectangle &rasterBBox, FeatureStats &stats );
 
     //! Tests whether a pixel's value should be included in the result
-    bool validPixel( float value ) const;
+    bool validPixel( double value ) const;
 
     QString getUniqueFieldName( const QString &fieldName, const QList<QgsField> &newFields );
 
@@ -153,8 +153,6 @@ class ANALYSIS_EXPORT QgsZonalStatistics
     int mRasterBand = 0;
     QgsVectorLayer *mPolygonLayer = nullptr;
     QString mAttributePrefix;
-    //! The nodata value of the input layer
-    float mInputNodataValue = -1;
     Statistics mStatistics = QgsZonalStatistics::All;
 };
 

tests/src/analysis/testqgszonalstatistics.cpp

@@ -40,10 +40,12 @@ class TestQgsZonalStatistics : public QObject
 
     void testStatistics();
     void testReprojection();
+    void testNoData();
 
   private:
     QgsVectorLayer *mVectorLayer = nullptr;
     QgsRasterLayer *mRasterLayer = nullptr;
+    QString mTempPath;
 };
 
 void TestQgsZonalStatistics::initTestCase()
@@ -54,19 +56,19 @@ void TestQgsZonalStatistics::initTestCase()
 
   QString myDataPath( TEST_DATA_DIR ); //defined in CmakeLists.txt
   QString myTestDataPath = myDataPath + "/zonalstatistics/";
-  QString myTempPath = QDir::tempPath() + '/';
+  mTempPath = QDir::tempPath() + '/';
 
   // copy test data to temp directory
   QDir testDir( myTestDataPath );
   QStringList files = testDir.entryList( QDir::Files | QDir::NoDotAndDotDot );
   for ( int i = 0; i < files.size(); ++i )
   {
-    QFile::remove( myTempPath + files.at( i ) );
-    QVERIFY( QFile::copy( myTestDataPath + files.at( i ), myTempPath + files.at( i ) ) );
+    QFile::remove( mTempPath + files.at( i ) );
+    QVERIFY( QFile::copy( myTestDataPath + files.at( i ), mTempPath + files.at( i ) ) );
   }
 
-  mVectorLayer = new QgsVectorLayer( myTempPath + "polys.shp", QStringLiteral( "poly" ), QStringLiteral( "ogr" ) );
-  mRasterLayer = new QgsRasterLayer( myTempPath + "edge_problem.asc", QStringLiteral( "raster" ), QStringLiteral( "gdal" ) );
+  mVectorLayer = new QgsVectorLayer( mTempPath + "polys.shp", QStringLiteral( "poly" ), QStringLiteral( "ogr" ) );
+  mRasterLayer = new QgsRasterLayer( mTempPath + "edge_problem.asc", QStringLiteral( "raster" ), QStringLiteral( "gdal" ) );
   QgsProject::instance()->addMapLayers(
     QList<QgsMapLayer *>() << mVectorLayer << mRasterLayer );
 }
@@ -246,5 +248,59 @@ void TestQgsZonalStatistics::testReprojection()
   QCOMPARE( f.attribute( "variance" ).toDouble(), 0.13888888888889 );
 }
 
+void TestQgsZonalStatistics::testNoData()
+{
+  QString myDataPath( TEST_DATA_DIR ); //defined in CmakeLists.txt
+  QString myTestDataPath = myDataPath + "/zonalstatistics/";
+
+  // test that zonal stats respects no data and user set no data values
+  std::unique_ptr< QgsRasterLayer > rasterLayer = qgis::make_unique< QgsRasterLayer >( myTestDataPath + "raster.tif", QStringLiteral( "raster" ), QStringLiteral( "gdal" ) );
+  std::unique_ptr< QgsVectorLayer > vectorLayer = qgis::make_unique< QgsVectorLayer >( mTempPath + "polys2.shp", QStringLiteral( "poly" ), QStringLiteral( "ogr" ) );
+
+  QgsZonalStatistics zs( vectorLayer.get(), rasterLayer.get(), QStringLiteral( "n" ), 1, QgsZonalStatistics::All );
+  zs.calculateStatistics( nullptr );
+
+  QgsFeature f;
+  QgsFeatureRequest request;
+  QgsFeatureIterator it = vectorLayer->getFeatures( request );
+  bool fetched = it.nextFeature( f );
+  QVERIFY( fetched );
+  QCOMPARE( f.attribute( "ncount" ).toDouble(), 16.0 );
+  QCOMPARE( f.attribute( "nsum" ).toDouble(), 13428.0 );
+
+  fetched = it.nextFeature( f );
+  QVERIFY( fetched );
+  QCOMPARE( f.attribute( "ncount" ).toDouble(), 103.0 );
+  QCOMPARE( f.attribute( "nsum" ).toDouble(), 90536.0 );
+
+  fetched = it.nextFeature( f );
+  QVERIFY( fetched );
+  QCOMPARE( f.attribute( "ncount" ).toDouble(), 0.0 );
+  QCOMPARE( f.attribute( "nsum" ).toDouble(), 0.0 );
+
+  // with user no data
+  rasterLayer->dataProvider()->setUserNoDataValue( 1, QgsRasterRangeList() << QgsRasterRange( 842, 852 )
+      << QgsRasterRange( 877, 891 ) );
+
+  zs = QgsZonalStatistics( vectorLayer.get(), rasterLayer.get(), QStringLiteral( "un" ), 1, QgsZonalStatistics::All );
+  zs.calculateStatistics( nullptr );
+
+  it = vectorLayer->getFeatures( request );
+  fetched = it.nextFeature( f );
+  QVERIFY( fetched );
+  QCOMPARE( f.attribute( "uncount" ).toDouble(), 8.0 );
+  QCOMPARE( f.attribute( "unsum" ).toDouble(), 6652.0 );
+
+  fetched = it.nextFeature( f );
+  QVERIFY( fetched );
+  QCOMPARE( f.attribute( "uncount" ).toDouble(), 52.0 );
+  QCOMPARE( f.attribute( "unsum" ).toDouble(), 45374.0 );
+
+  fetched = it.nextFeature( f );
+  QVERIFY( fetched );
+  QCOMPARE( f.attribute( "uncount" ).toDouble(), 0.0 );
+  QCOMPARE( f.attribute( "unsum" ).toDouble(), 0.0 );
+}
+
 QGSTEST_MAIN( TestQgsZonalStatistics )
 #include "testqgszonalstatistics.moc"

tests/testdata/zonalstatistics/edge_problem.asc.aux.xml

@@ -17,6 +17,14 @@
         <Approximate>0</Approximate>
         <HistCounts>4|0|0|0|0|0|0|0|0|0|0|8</HistCounts>
       </HistItem>
+      <HistItem>
+        <HistMin>-0.25</HistMin>
+        <HistMax>1.25</HistMax>
+        <BucketCount>2</BucketCount>
+        <IncludeOutOfRange>0</IncludeOutOfRange>
+        <Approximate>0</Approximate>
+        <HistCounts>4|8</HistCounts>
+      </HistItem>
     </Histograms>
     <Metadata>
       <MDI key="STATISTICS_MAXIMUM">1</MDI>

tests/testdata/zonalstatistics/polys2.cpg

@@ -0,0 +1 @@
+UTF-8

tests/testdata/zonalstatistics/polys2.dbf

@@ -0,0 +1 @@
+PROJCS["ED50_UTM_zone_30N",GEOGCS["GCS_European_1950",DATUM["D_European_1950",SPHEROID["International_1924",6378388,297]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-3],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]

tests/testdata/zonalstatistics/polys2.prj

@@ -0,0 +1 @@
+PROJCS["ED50 / UTM zone 30N",GEOGCS["ED50",DATUM["European_Datum_1950",SPHEROID["International 1924",6378388,297,AUTHORITY["EPSG","7022"]],TOWGS84[-87,-98,-121,0,0,0,0],AUTHORITY["EPSG","6230"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4230"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-3],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","23030"]]

tests/testdata/zonalstatistics/polys2.qpj

@@ -0,0 +1,10 @@
+<PAMDataset>
+  <PAMRasterBand band="1">
+    <Metadata>
+      <MDI key="STATISTICS_MAXIMUM">899</MDI>
+      <MDI key="STATISTICS_MEAN">865.86666666667</MDI>
+      <MDI key="STATISTICS_MINIMUM">826</MDI>
+      <MDI key="STATISTICS_STDDEV">17.808206597584</MDI>
+    </Metadata>
+  </PAMRasterBand>
+</PAMDataset>