Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Fix calculation of zonal stats when source contains nodata or nan
pixels (fix #11135)
  • Loading branch information
nyalldawson committed Apr 16, 2015
1 parent 6a46f46 commit 213a22b
Show file tree
Hide file tree
Showing 2 changed files with 24 additions and 11 deletions.
31 changes: 20 additions & 11 deletions src/analysis/vector/qgszonalstatistics.cpp
Expand Up @@ -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 );
}
Expand All @@ -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 )
{
Expand All @@ -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();
Expand Down
4 changes: 4 additions & 0 deletions src/analysis/vector/qgszonalstatistics.h
Expand Up @@ -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;
Expand All @@ -63,6 +66,7 @@ class ANALYSIS_EXPORT QgsZonalStatistics
QString mAttributePrefix;
/**The nodata value of the input layer*/
float mInputNodataValue;

};

#endif // QGSZONALSTATISTICS_H

0 comments on commit 213a22b

Please sign in to comment.