Fix zonal statistics calculations when pixel size is greater than

polygon size Fixes #17159
qgis · Jun 8, 2018 · 83d44b9 · 83d44b9
1 parent 985aee6
commit 83d44b9
Show file tree

Hide file tree

Showing 8 changed files with 58 additions and 16 deletions.
diff --git a/src/analysis/vector/qgszonalstatistics.cpp b/src/analysis/vector/qgszonalstatistics.cpp
@@ -456,25 +456,23 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( const QgsGeometry &p
     for ( int j = 0; j < nCellsX; ++j )
     {
       double pixelValue = block->value( i, j );
-      if ( !validPixel( pixelValue ) || block->isNoData( i, j ) )
-      {
-        continue;
-      }
-
-      pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) );
-      // GEOS intersects tests on prepared geometry is MAGNITUDES faster than calculating the intersection itself,
-      // so we first test to see if there IS an intersection before doing the actual calculation
-      if ( !pixelRectGeometry.isNull() && polyEngine->intersects( pixelRectGeometry.constGet() ) )
+      if ( validPixel( pixelValue ) && !block->isNoData( i, j ) )
       {
-        //intersection
-        QgsGeometry intersectGeometry = pixelRectGeometry.intersection( poly );
-        if ( !intersectGeometry.isEmpty() )
+        pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) );
+        // GEOS intersects tests on prepared geometry is MAGNITUDES faster than calculating the intersection itself,
+        // so we first test to see if there IS an intersection before doing the actual calculation
+        if ( !pixelRectGeometry.isNull() && polyEngine->intersects( pixelRectGeometry.constGet() ) )
         {
-          double intersectionArea = intersectGeometry.area();
-          if ( intersectionArea >= 0.0 )
+          //intersection
+          QgsGeometry intersectGeometry = pixelRectGeometry.intersection( poly );
+          if ( !intersectGeometry.isEmpty() )
           {
-            weight = intersectionArea / pixelArea;
-            stats.addValue( pixelValue, weight );
+            double intersectionArea = intersectGeometry.area();
+            if ( intersectionArea >= 0.0 )
+            {
+              weight = intersectionArea / pixelArea;
+              stats.addValue( pixelValue, weight );
+            }
           }
         }
       }

diff --git a/tests/src/analysis/testqgszonalstatistics.cpp b/tests/src/analysis/testqgszonalstatistics.cpp
@@ -41,6 +41,7 @@ class TestQgsZonalStatistics : public QObject
     void testStatistics();
     void testReprojection();
     void testNoData();
+    void testSmallPolygons();
 
   private:
     QgsVectorLayer *mVectorLayer = nullptr;
@@ -302,5 +303,45 @@ void TestQgsZonalStatistics::testNoData()
   QCOMPARE( f.attribute( "unsum" ).toDouble(), 0.0 );
 }
 
+void TestQgsZonalStatistics::testSmallPolygons()
+{
+  QString myDataPath( TEST_DATA_DIR ); //defined in CmakeLists.txt
+  QString myTestDataPath = myDataPath + "/zonalstatistics/";
+
+  // test that zonal stats works ok with polygons much smaller than pixel size
+  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 + "small_polys.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 );
+  QGSCOMPARENEAR( f.attribute( "ncount" ).toDouble(), 0.698248, 0.001 );
+  QGSCOMPARENEAR( f.attribute( "nsum" ).toDouble(), 588.711, 0.001 );
+  QCOMPARE( f.attribute( "nmin" ).toDouble(), 826.0 );
+  QCOMPARE( f.attribute( "nmax" ).toDouble(), 851.0 );
+  QGSCOMPARENEAR( f.attribute( "nmean" ).toDouble(), 843.125292, 0.001 );
+
+  fetched = it.nextFeature( f );
+  QVERIFY( fetched );
+  QGSCOMPARENEAR( f.attribute( "ncount" ).toDouble(), 0.240808, 0.001 );
+  QGSCOMPARENEAR( f.attribute( "nsum" ).toDouble(), 208.030921, 0.001 );
+  QCOMPARE( f.attribute( "nmin" ).toDouble(), 859.0 );
+  QCOMPARE( f.attribute( "nmax" ).toDouble(), 868.0 );
+  QGSCOMPARENEAR( f.attribute( "nmean" ).toDouble(), 863.887500, 0.001 );
+
+  fetched = it.nextFeature( f );
+  QVERIFY( fetched );
+  QGSCOMPARENEAR( f.attribute( "ncount" ).toDouble(), 0.259522, 0.001 );
+  QGSCOMPARENEAR( f.attribute( "nsum" ).toDouble(), 224.300747, 0.001 );
+  QCOMPARE( f.attribute( "nmin" ).toDouble(), 851.0 );
+  QCOMPARE( f.attribute( "nmax" ).toDouble(), 872.0 );
+  QGSCOMPARENEAR( f.attribute( "nmean" ).toDouble(), 864.285638, 0.001 );
+}
+
 QGSTEST_MAIN( TestQgsZonalStatistics )
 #include "testqgszonalstatistics.moc"
diff --git a/tests/testdata/zonalstatistics/small_polys.cpg b/tests/testdata/zonalstatistics/small_polys.cpg
@@ -0,0 +1 @@
+UTF-8
diff --git a/tests/testdata/zonalstatistics/small_polys.dbf b/tests/testdata/zonalstatistics/small_polys.dbf
diff --git a/tests/testdata/zonalstatistics/small_polys.prj b/tests/testdata/zonalstatistics/small_polys.prj
@@ -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]]
diff --git a/tests/testdata/zonalstatistics/small_polys.qpj b/tests/testdata/zonalstatistics/small_polys.qpj
@@ -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"]]
diff --git a/tests/testdata/zonalstatistics/small_polys.shp b/tests/testdata/zonalstatistics/small_polys.shp
diff --git a/tests/testdata/zonalstatistics/small_polys.shx b/tests/testdata/zonalstatistics/small_polys.shx