Skip to content

Commit 83d44b9

Browse files
nyalldawsonnyalldawson
committedJun 8, 2018
Fix zonal statistics calculations when pixel size is greater than
polygon size Fixes #17159
·
ltr-3_40final-3_2_0
1 parent 985aee6 commit 83d44b9

File tree

8 files changed

+58
-16
lines changed

8 files changed

+58
-16
lines changed
 

‎src/analysis/vector/qgszonalstatistics.cpp

Lines changed: 14 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -456,25 +456,23 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( const QgsGeometry &p
456456
for ( int j = 0; j < nCellsX; ++j )
457457
{
458458
double pixelValue = block->value( i, j );
459-
if ( !validPixel( pixelValue ) || block->isNoData( i, j ) )
460-
{
461-
continue;
462-
}
463-
464-
pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) );
465-
// GEOS intersects tests on prepared geometry is MAGNITUDES faster than calculating the intersection itself,
466-
// so we first test to see if there IS an intersection before doing the actual calculation
467-
if ( !pixelRectGeometry.isNull() && polyEngine->intersects( pixelRectGeometry.constGet() ) )
459+
if ( validPixel( pixelValue ) && !block->isNoData( i, j ) )
468460
{
469-
//intersection
470-
QgsGeometry intersectGeometry = pixelRectGeometry.intersection( poly );
471-
if ( !intersectGeometry.isEmpty() )
461+
pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) );
462+
// GEOS intersects tests on prepared geometry is MAGNITUDES faster than calculating the intersection itself,
463+
// so we first test to see if there IS an intersection before doing the actual calculation
464+
if ( !pixelRectGeometry.isNull() && polyEngine->intersects( pixelRectGeometry.constGet() ) )
472465
{
473-
double intersectionArea = intersectGeometry.area();
474-
if ( intersectionArea >= 0.0 )
466+
//intersection
467+
QgsGeometry intersectGeometry = pixelRectGeometry.intersection( poly );
468+
if ( !intersectGeometry.isEmpty() )
475469
{
476-
weight = intersectionArea / pixelArea;
477-
stats.addValue( pixelValue, weight );
470+
double intersectionArea = intersectGeometry.area();
471+
if ( intersectionArea >= 0.0 )
472+
{
473+
weight = intersectionArea / pixelArea;
474+
stats.addValue( pixelValue, weight );
475+
}
478476
}
479477
}
480478
}

‎tests/src/analysis/testqgszonalstatistics.cpp

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,7 @@ class TestQgsZonalStatistics : public QObject
4141
void testStatistics();
4242
void testReprojection();
4343
void testNoData();
44+
void testSmallPolygons();
4445

4546
private:
4647
QgsVectorLayer *mVectorLayer = nullptr;
@@ -302,5 +303,45 @@ void TestQgsZonalStatistics::testNoData()
302303
QCOMPARE( f.attribute( "unsum" ).toDouble(), 0.0 );
303304
}
304305

306+
void TestQgsZonalStatistics::testSmallPolygons()
307+
{
308+
QString myDataPath( TEST_DATA_DIR ); //defined in CmakeLists.txt
309+
QString myTestDataPath = myDataPath + "/zonalstatistics/";
310+
311+
// test that zonal stats works ok with polygons much smaller than pixel size
312+
std::unique_ptr< QgsRasterLayer > rasterLayer = qgis::make_unique< QgsRasterLayer >( myTestDataPath + "raster.tif", QStringLiteral( "raster" ), QStringLiteral( "gdal" ) );
313+
std::unique_ptr< QgsVectorLayer > vectorLayer = qgis::make_unique< QgsVectorLayer >( mTempPath + "small_polys.shp", QStringLiteral( "poly" ), QStringLiteral( "ogr" ) );
314+
315+
QgsZonalStatistics zs( vectorLayer.get(), rasterLayer.get(), QStringLiteral( "n" ), 1, QgsZonalStatistics::All );
316+
zs.calculateStatistics( nullptr );
317+
318+
QgsFeature f;
319+
QgsFeatureRequest request;
320+
QgsFeatureIterator it = vectorLayer->getFeatures( request );
321+
bool fetched = it.nextFeature( f );
322+
QVERIFY( fetched );
323+
QGSCOMPARENEAR( f.attribute( "ncount" ).toDouble(), 0.698248, 0.001 );
324+
QGSCOMPARENEAR( f.attribute( "nsum" ).toDouble(), 588.711, 0.001 );
325+
QCOMPARE( f.attribute( "nmin" ).toDouble(), 826.0 );
326+
QCOMPARE( f.attribute( "nmax" ).toDouble(), 851.0 );
327+
QGSCOMPARENEAR( f.attribute( "nmean" ).toDouble(), 843.125292, 0.001 );
328+
329+
fetched = it.nextFeature( f );
330+
QVERIFY( fetched );
331+
QGSCOMPARENEAR( f.attribute( "ncount" ).toDouble(), 0.240808, 0.001 );
332+
QGSCOMPARENEAR( f.attribute( "nsum" ).toDouble(), 208.030921, 0.001 );
333+
QCOMPARE( f.attribute( "nmin" ).toDouble(), 859.0 );
334+
QCOMPARE( f.attribute( "nmax" ).toDouble(), 868.0 );
335+
QGSCOMPARENEAR( f.attribute( "nmean" ).toDouble(), 863.887500, 0.001 );
336+
337+
fetched = it.nextFeature( f );
338+
QVERIFY( fetched );
339+
QGSCOMPARENEAR( f.attribute( "ncount" ).toDouble(), 0.259522, 0.001 );
340+
QGSCOMPARENEAR( f.attribute( "nsum" ).toDouble(), 224.300747, 0.001 );
341+
QCOMPARE( f.attribute( "nmin" ).toDouble(), 851.0 );
342+
QCOMPARE( f.attribute( "nmax" ).toDouble(), 872.0 );
343+
QGSCOMPARENEAR( f.attribute( "nmean" ).toDouble(), 864.285638, 0.001 );
344+
}
345+
305346
QGSTEST_MAIN( TestQgsZonalStatistics )
306347
#include "testqgszonalstatistics.moc"

‎tests/testdata/zonalstatistics/small_polys.cpg

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
UTF-8

‎tests/testdata/zonalstatistics/small_polys.dbf

99 Bytes
Binary file not shown.

‎tests/testdata/zonalstatistics/small_polys.prj

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
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/small_polys.qpj

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
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/small_polys.shp

508 Bytes
Binary file not shown.

‎tests/testdata/zonalstatistics/small_polys.shx

124 Bytes
Binary file not shown.

0 commit comments

Comments
 (0)
Please sign in to comment.