Skip to content
Permalink
Browse files

Fix zonal statistics calculations when pixel size is greater than

polygon size

Fixes #17159
  • Loading branch information
nyalldawson committed Jun 8, 2018
1 parent 985aee6 commit 83d44b9fe975dbfcd1007b347f02712c08d0baff
@@ -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 );
}
}
}
}
@@ -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"
@@ -0,0 +1 @@
UTF-8
Binary file not shown.
@@ -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]]
@@ -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"]]
Binary file not shown.
Binary file not shown.

0 comments on commit 83d44b9

Please sign in to comment.
You can’t perform that action at this time.