Skip to content

Commit

Permalink
use QgsRasterBlock instead of GDAL in zonal statistics (addresses #8736)
Browse files Browse the repository at this point in the history 
This should make zonal statistics usable rasters which come from
other providers, e.g. WCS.
  • Loading branch information
@alexbruy
alexbruy committed Jan 28, 2017
1 parent 290758a commit 278913b
Show file tree
Hide file tree
Showing 6 changed files with 61 additions and 74 deletions.
2 changes: 1 addition & 1 deletion python/analysis/vector/qgszonalstatistics.sip
View file
Expand Up @@ -30,7 +30,7 @@ class QgsZonalStatistics

typedef QFlags<QgsZonalStatistics::Statistic> Statistics;

QgsZonalStatistics( QgsVectorLayer* polygonLayer, const QString& rasterFile, const QString& attributePrefix = "", int rasterBand = 1,
QgsZonalStatistics( QgsVectorLayer* polygonLayer, QgsRasterLayer* rasterLayer, const QString& attributePrefix = "", int rasterBand = 1,
QgsZonalStatistics::Statistics stats = QgsZonalStatistics::Statistics( QgsZonalStatistics::Count | QgsZonalStatistics::Sum | QgsZonalStatistics::Mean) );

/** Starts the calculation
Expand Down
3 changes: 2 additions & 1 deletion python/plugins/processing/algs/qgis/ZonalStatisticsQgis.py
View file
Expand Up @@ -100,14 +100,15 @@ def processAlgorithm(self, feedback):
st = self.getParameterValue(self.STATISTICS)

vectorLayer = dataobjects.getObjectFromUri(vectorPath)
rasterLayer = dataobjects.getObjectFromUri(rasterPath)

keys = list(self.STATS.keys())
selectedStats = 0
for i in st:
selectedStats |= self.STATS[keys[i]]

zs = QgsZonalStatistics(vectorLayer,
rasterPath,
rasterLayer,
columnPrefix,
bandNumber,
selectedStats)
Expand Down
101 changes: 42 additions & 59 deletions src/analysis/vector/qgszonalstatistics.cpp
View file
Expand Up @@ -20,16 +20,17 @@
#include "qgsgeometry.h"
#include "qgsvectordataprovider.h"
#include "qgsvectorlayer.h"
#include "qgsrasterdataprovider.h"
#include "qgsrasterlayer.h"
#include "qgsrasterblock.h"
#include "qmath.h"
#include "gdal.h"
#include "cpl_string.h"
#include "qgslogger.h"

#include <QProgressDialog>
#include <QFile>

QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer* polygonLayer, const QString& rasterFile, const QString& attributePrefix, int rasterBand, Statistics stats )
: mRasterFilePath( rasterFile )
QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer* polygonLayer, QgsRasterLayer* rasterLayer, const QString& attributePrefix, int rasterBand, Statistics stats )
: mRasterLayer( rasterLayer )
, mRasterBand( rasterBand )
, mPolygonLayer( polygonLayer )
, mAttributePrefix( attributePrefix )
Expand All @@ -40,7 +41,8 @@ QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer* polygonLayer, const QStr
}

QgsZonalStatistics::QgsZonalStatistics()
: mRasterBand( 0 )
: mRasterLayer( nullptr )
, mRasterBand( 0 )
, mPolygonLayer( nullptr )
, mInputNodataValue( -1 )
, mStatistics( QgsZonalStatistics::All )
Expand All @@ -61,49 +63,33 @@ int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
return 2;
}

//open the raster layer and the raster band
GDALAllRegister();
GDALDatasetH inputDataset = GDALOpen( mRasterFilePath.toUtf8().constData(), GA_ReadOnly );
if ( !inputDataset )
if ( !mRasterLayer )
{
return 3;
}

if ( GDALGetRasterCount( inputDataset ) < ( mRasterBand - 1 ) )
if ( mRasterLayer->bandCount() < mRasterBand )
{
GDALClose( inputDataset );
return 4;
}

GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset, mRasterBand );
if ( !rasterBand )
{
GDALClose( inputDataset );
return 5;
}
mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, nullptr );
mRasterProvider = mRasterLayer->dataProvider();
mInputNodataValue = mRasterProvider->sourceNoDataValue( mRasterBand );

//get geometry info about raster layer
int nCellsXGDAL = GDALGetRasterXSize( inputDataset );
int nCellsYGDAL = GDALGetRasterYSize( inputDataset );
double geoTransform[6];
if ( GDALGetGeoTransform( inputDataset, geoTransform ) != CE_None )
{
GDALClose( inputDataset );
return 6;
}
double cellsizeX = geoTransform[1];
int nCellsXProvider = mRasterProvider->xSize();
int nCellsYProvider = mRasterProvider->ySize();
double cellsizeX = mRasterLayer->rasterUnitsPerPixelX();
if ( cellsizeX < 0 )
{
cellsizeX = -cellsizeX;
}
double cellsizeY = geoTransform[5];
double cellsizeY = mRasterLayer->rasterUnitsPerPixelY();
if ( cellsizeY < 0 )
{
cellsizeY = -cellsizeY;
}
QgsRectangle rasterBBox( geoTransform[0], geoTransform[3] - ( nCellsYGDAL * cellsizeY ),
geoTransform[0] + ( nCellsXGDAL * cellsizeX ), geoTransform[3] );
QgsRectangle rasterBBox = mRasterProvider->extent();

//add the new fields to the provider
QList<QgsField> newFieldList;
Expand Down Expand Up @@ -223,7 +209,6 @@ int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
p->setMaximum( featureCount );
}


//iterate over each polygon
QgsFeatureRequest request;
request.setSubsetOfAttributes( QgsAttributeList() );
Expand Down Expand Up @@ -273,22 +258,22 @@ int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
}

//avoid access to cells outside of the raster (may occur because of rounding)
if (( offsetX + nCellsX ) > nCellsXGDAL )
if (( offsetX + nCellsX ) > nCellsXProvider )
{
nCellsX = nCellsXGDAL - offsetX;
nCellsX = nCellsXProvider - offsetX;
}
if (( offsetY + nCellsY ) > nCellsYGDAL )
if (( offsetY + nCellsY ) > nCellsYProvider )
{
nCellsY = nCellsYGDAL - offsetY;
nCellsY = nCellsYProvider - offsetY;
}

statisticsFromMiddlePointTest( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
statisticsFromMiddlePointTest( featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
rasterBBox, featureStats );

if ( featureStats.count <= 1 )
{
//the cell resolution is probably larger than the polygon area. We switch to precise pixel - polygon intersection in this case
statisticsFromPreciseIntersection( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
statisticsFromPreciseIntersection( featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
rasterBBox, featureStats );
}

Expand Down Expand Up @@ -366,7 +351,6 @@ int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
p->setValue( featureCount );
}

GDALClose( inputDataset );
mPolygonLayer->updateFields();

if ( p && p->wasCanceled() )
Expand Down Expand Up @@ -404,12 +388,11 @@ int QgsZonalStatistics::cellInfoForBBox( const QgsRectangle& rasterBBox, const Q
return 0;
}

void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, const QgsGeometry& poly, int pixelOffsetX,
void QgsZonalStatistics::statisticsFromMiddlePointTest( const QgsGeometry& poly, int pixelOffsetX,
int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, FeatureStats &stats )
{
double cellCenterX, cellCenterY;

float* scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );
cellCenterY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
stats.reset();

Expand All @@ -430,17 +413,16 @@ void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, const QgsGeo
GEOSCoordSequence* cellCenterCoords = nullptr;
GEOSGeometry* currentCellCenter = nullptr;

QgsRectangle featureBBox = poly.boundingBox().intersect( &rasterBBox );
QgsRectangle intersectBBox = rasterBBox.intersect( &featureBBox );

QgsRasterBlock* block = mRasterProvider->block( mRasterBand, intersectBBox, nCellsX, nCellsY );
for ( int i = 0; i < nCellsY; ++i )
{
if ( GDALRasterIO( band, GF_Read, pixelOffsetX, pixelOffsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 )
!= CPLE_None )
{
continue;
}
cellCenterX = rasterBBox.xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;
for ( int j = 0; j < nCellsX; ++j )
{
if ( validPixel( scanLine[j] ) )
if ( validPixel( block->value( i, j ) ) )
{
GEOSGeom_destroy_r( geosctxt, currentCellCenter );
cellCenterCoords = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
Expand All @@ -449,45 +431,46 @@ void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, const QgsGeo
currentCellCenter = GEOSGeom_createPoint_r( geosctxt, cellCenterCoords );
if ( GEOSPreparedContains_r( geosctxt, polyGeosPrepared, currentCellCenter ) )
{
stats.addValue( scanLine[j] );
stats.addValue( block->value( i, j ) );
}
}
cellCenterX += cellSizeX;
}
cellCenterY -= cellSizeY;
}

GEOSGeom_destroy_r( geosctxt, currentCellCenter );
CPLFree( scanLine );
GEOSPreparedGeom_destroy_r( geosctxt, polyGeosPrepared );
GEOSGeom_destroy_r( geosctxt, polyGeos );
delete block;
}

void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, const QgsGeometry& poly, int pixelOffsetX,
void QgsZonalStatistics::statisticsFromPreciseIntersection( const QgsGeometry& poly, int pixelOffsetX,
int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, FeatureStats &stats )
{
stats.reset();

double currentY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
float* pixelData = ( float * ) CPLMalloc( sizeof( float ) );
QgsGeometry pixelRectGeometry;

double hCellSizeX = cellSizeX / 2.0;
double hCellSizeY = cellSizeY / 2.0;
double pixelArea = cellSizeX * cellSizeY;
double weight = 0;

for ( int row = 0; row < nCellsY; ++row )
QgsRectangle featureBBox = poly.boundingBox().intersect( &rasterBBox );
QgsRectangle intersectBBox = rasterBBox.intersect( &featureBBox );

QgsRasterBlock* block = mRasterProvider->block( mRasterBand, intersectBBox, nCellsX, nCellsY );
for ( int i = 0; i < nCellsY; ++i )
{
double currentX = rasterBBox.xMinimum() + cellSizeX / 2.0 + pixelOffsetX * cellSizeX;
for ( int col = 0; col < nCellsX; ++col )
for ( int j = 0; j < nCellsX; ++j )
{
if ( GDALRasterIO( band, GF_Read, pixelOffsetX + col, pixelOffsetY + row, nCellsX, 1, pixelData, 1, 1, GDT_Float32, 0, 0 ) != CE_None )
if ( !validPixel( block->value( i, j ) ) )
{
QgsDebugMsg( "Raster IO Error" );
}

if ( !validPixel( *pixelData ) )
continue;
}

pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) );
if ( !pixelRectGeometry.isEmpty() )
Expand All @@ -500,7 +483,7 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, const Qg
if ( intersectionArea >= 0.0 )
{
weight = intersectionArea / pixelArea;
stats.addValue( *pixelData, weight );
stats.addValue( block->value( i, j ), weight );
}
}
pixelRectGeometry = QgsGeometry();
Expand All @@ -509,7 +492,7 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, const Qg
}
currentY -= cellSizeY;
}
CPLFree( pixelData );
delete block;
}

bool QgsZonalStatistics::validPixel( float value ) const
Expand Down
11 changes: 7 additions & 4 deletions src/analysis/vector/qgszonalstatistics.h
View file
Expand Up @@ -26,6 +26,8 @@

class QgsGeometry;
class QgsVectorLayer;
class QgsRasterLayer;
class QgsRasterDataProvider;
class QProgressDialog;
class QgsRectangle;
class QgsField;
Expand Down Expand Up @@ -57,7 +59,7 @@ class ANALYSIS_EXPORT QgsZonalStatistics
/**
* Constructor for QgsZonalStatistics.
*/
QgsZonalStatistics( QgsVectorLayer* polygonLayer, const QString& rasterFile, const QString& attributePrefix = "", int rasterBand = 1,
QgsZonalStatistics( QgsVectorLayer* polygonLayer, QgsRasterLayer* rasterLayer, const QString& attributePrefix = "", int rasterBand = 1,
Statistics stats = Statistics( Count | Sum | Mean ) );

/** Starts the calculation
Expand Down Expand Up @@ -114,19 +116,20 @@ class ANALYSIS_EXPORT QgsZonalStatistics
int& offsetX, int& offsetY, int& nCellsX, int& nCellsY ) const;

//! Returns statistics by considering the pixels where the center point is within the polygon (fast)
void statisticsFromMiddlePointTest( void* band, const QgsGeometry& poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY,
void statisticsFromMiddlePointTest( const QgsGeometry& poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY,
double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, FeatureStats& stats );

//! Returns statistics with precise pixel - polygon intersection test (slow)
void statisticsFromPreciseIntersection( void* band, const QgsGeometry& poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY,
void statisticsFromPreciseIntersection( const QgsGeometry& poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY,
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;

QString getUniqueFieldName( const QString& fieldName, const QList<QgsField>& newFields );

QString mRasterFilePath;
QgsRasterLayer* mRasterLayer;
QgsRasterDataProvider* mRasterProvider;
//! Raster band to calculate statistics from (defaults to 1)
int mRasterBand;
QgsVectorLayer* mPolygonLayer;
Expand Down
12 changes: 6 additions & 6 deletions tests/src/analysis/testqgszonalstatistics.cpp
View file
Expand Up @@ -19,6 +19,7 @@
#include "qgsapplication.h"
#include "qgsfeatureiterator.h"
#include "qgsvectorlayer.h"
#include "qgsrasterlayer.h"
#include "qgszonalstatistics.h"
#include "qgsproject.h"

Expand All @@ -42,7 +43,7 @@ class TestQgsZonalStatistics : public QObject

private:
QgsVectorLayer* mVectorLayer;
QString mRasterPath;
QgsRasterLayer* mRasterLayer;
};

TestQgsZonalStatistics::TestQgsZonalStatistics()
Expand Down Expand Up @@ -71,10 +72,9 @@ void TestQgsZonalStatistics::initTestCase()
}

mVectorLayer = new QgsVectorLayer( myTempPath + "polys.shp", QStringLiteral( "poly" ), QStringLiteral( "ogr" ) );
mRasterLayer = new QgsRasterLayer( myTempPath + "edge_problem.asc", QStringLiteral( "raster" ), QStringLiteral( "gdal" ) );
QgsProject::instance()->addMapLayers(
QList<QgsMapLayer *>() << mVectorLayer );

mRasterPath = myTempPath + "edge_problem.asc";
QList<QgsMapLayer *>() << mVectorLayer << mRasterLayer );
}

void TestQgsZonalStatistics::cleanupTestCase()
Expand All @@ -84,7 +84,7 @@ void TestQgsZonalStatistics::cleanupTestCase()

void TestQgsZonalStatistics::testStatistics()
{
QgsZonalStatistics zs( mVectorLayer, mRasterPath, QLatin1String( "" ), 1, QgsZonalStatistics::All );
QgsZonalStatistics zs( mVectorLayer, mRasterLayer, QLatin1String( "" ), 1, QgsZonalStatistics::All );
zs.calculateStatistics( nullptr );

QgsFeature f;
Expand Down Expand Up @@ -135,7 +135,7 @@ void TestQgsZonalStatistics::testStatistics()
QCOMPARE( f.attribute( "variety" ).toDouble(), 2.0 );

// same with long prefix to ensure that field name truncation handled correctly
QgsZonalStatistics zsl( mVectorLayer, mRasterPath, QStringLiteral( "myqgis2_" ), 1, QgsZonalStatistics::All );
QgsZonalStatistics zsl( mVectorLayer, mRasterLayer, QStringLiteral( "myqgis2_" ), 1, QgsZonalStatistics::All );
zsl.calculateStatistics( nullptr );

request.setFilterFid( 0 );
Expand Down
6 changes: 3 additions & 3 deletions tests/src/python/test_qgszonalstatistics.py
View file
Expand Up @@ -15,7 +15,7 @@
import qgis # NOQA

from qgis.PyQt.QtCore import QDir, QFile
from qgis.core import QgsVectorLayer, QgsFeature, QgsFeatureRequest
from qgis.core import QgsVectorLayer, QgsRasterLayer, QgsFeature, QgsFeatureRequest
from qgis.analysis import QgsZonalStatistics

from qgis.testing import start_app, unittest
Expand All @@ -38,8 +38,8 @@ def testStatistics(self):
QFile.copy(TEST_DATA_DIR + f, myTempPath + f)

myVector = QgsVectorLayer(myTempPath + "polys.shp", "poly", "ogr")
myRasterPath = myTempPath + "edge_problem.asc"
zs = QgsZonalStatistics(myVector, myRasterPath, "", 1, QgsZonalStatistics.All)
myRaster = QgsRasterLayer(myTempPath + "edge_problem.asc", "raster", "gdal")
zs = QgsZonalStatistics(myVector, myRaster, "", 1, QgsZonalStatistics.All)
zs.calculateStatistics(None)

feat = QgsFeature()
Expand Down

0 comments on commit 278913b

Please sign in to comment.