Skip to content

Commit 95091d9

Browse files
authored
Merge pull request #4062 from alexbruy/raster-provider
Use native QGIS raster API instead of GDAL API in zonal statistics
2 parents 02b7a55 + d681e9c commit 95091d9

File tree

7 files changed

+64
-74
lines changed

7 files changed

+64
-74
lines changed

doc/api_break.dox

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2061,6 +2061,9 @@ Triangulation {#qgis_api_break_3_0_Triangulation}
20612061

20622062
- forcedCrossBehaviour enum and its setter and getter have been renamed to ForcedCrossBehavior. Its members have been renamed: SnappingType_VERTICE: SnappingTypeVertex, DELETE_FIRST: DeleteFirst, INSERT_VERTEX: InsertVertex. <!--#spellok-->
20632063

2064+
QgsZonalStatistics {#qgis_api_break_3_0_QgsZonalStatistics}
2065+
------------------
2066+
- QgsZonalStatistics() сonstructor now accepts pointer to the QgsRasterLayer instance instead of path to the raster file
20642067

20652068

20662069
QGIS 2.6 {#qgis_api_break_2_6}

python/analysis/vector/qgszonalstatistics.sip

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -30,7 +30,7 @@ class QgsZonalStatistics
3030

3131
typedef QFlags<QgsZonalStatistics::Statistic> Statistics;
3232

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

3636
/** Starts the calculation

python/plugins/processing/algs/qgis/ZonalStatisticsQgis.py

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -100,14 +100,15 @@ def processAlgorithm(self, feedback):
100100
st = self.getParameterValue(self.STATISTICS)
101101

102102
vectorLayer = dataobjects.getObjectFromUri(vectorPath)
103+
rasterLayer = dataobjects.getObjectFromUri(rasterPath)
103104

104105
keys = list(self.STATS.keys())
105106
selectedStats = 0
106107
for i in st:
107108
selectedStats |= self.STATS[keys[i]]
108109

109110
zs = QgsZonalStatistics(vectorLayer,
110-
rasterPath,
111+
rasterLayer,
111112
columnPrefix,
112113
bandNumber,
113114
selectedStats)

src/analysis/vector/qgszonalstatistics.cpp

Lines changed: 42 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -20,16 +20,17 @@
2020
#include "qgsgeometry.h"
2121
#include "qgsvectordataprovider.h"
2222
#include "qgsvectorlayer.h"
23+
#include "qgsrasterdataprovider.h"
24+
#include "qgsrasterlayer.h"
25+
#include "qgsrasterblock.h"
2326
#include "qmath.h"
24-
#include "gdal.h"
25-
#include "cpl_string.h"
2627
#include "qgslogger.h"
2728

2829
#include <QProgressDialog>
2930
#include <QFile>
3031

31-
QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer* polygonLayer, const QString& rasterFile, const QString& attributePrefix, int rasterBand, Statistics stats )
32-
: mRasterFilePath( rasterFile )
32+
QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer* polygonLayer, QgsRasterLayer* rasterLayer, const QString& attributePrefix, int rasterBand, Statistics stats )
33+
: mRasterLayer( rasterLayer )
3334
, mRasterBand( rasterBand )
3435
, mPolygonLayer( polygonLayer )
3536
, mAttributePrefix( attributePrefix )
@@ -40,7 +41,8 @@ QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer* polygonLayer, const QStr
4041
}
4142

4243
QgsZonalStatistics::QgsZonalStatistics()
43-
: mRasterBand( 0 )
44+
: mRasterLayer( nullptr )
45+
, mRasterBand( 0 )
4446
, mPolygonLayer( nullptr )
4547
, mInputNodataValue( -1 )
4648
, mStatistics( QgsZonalStatistics::All )
@@ -61,49 +63,33 @@ int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
6163
return 2;
6264
}
6365

64-
//open the raster layer and the raster band
65-
GDALAllRegister();
66-
GDALDatasetH inputDataset = GDALOpen( mRasterFilePath.toUtf8().constData(), GA_ReadOnly );
67-
if ( !inputDataset )
66+
if ( !mRasterLayer )
6867
{
6968
return 3;
7069
}
7170

72-
if ( GDALGetRasterCount( inputDataset ) < ( mRasterBand - 1 ) )
71+
if ( mRasterLayer->bandCount() < mRasterBand )
7372
{
74-
GDALClose( inputDataset );
7573
return 4;
7674
}
7775

78-
GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset, mRasterBand );
79-
if ( !rasterBand )
80-
{
81-
GDALClose( inputDataset );
82-
return 5;
83-
}
84-
mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, nullptr );
76+
mRasterProvider = mRasterLayer->dataProvider();
77+
mInputNodataValue = mRasterProvider->sourceNoDataValue( mRasterBand );
8578

8679
//get geometry info about raster layer
87-
int nCellsXGDAL = GDALGetRasterXSize( inputDataset );
88-
int nCellsYGDAL = GDALGetRasterYSize( inputDataset );
89-
double geoTransform[6];
90-
if ( GDALGetGeoTransform( inputDataset, geoTransform ) != CE_None )
91-
{
92-
GDALClose( inputDataset );
93-
return 6;
94-
}
95-
double cellsizeX = geoTransform[1];
80+
int nCellsXProvider = mRasterProvider->xSize();
81+
int nCellsYProvider = mRasterProvider->ySize();
82+
double cellsizeX = mRasterLayer->rasterUnitsPerPixelX();
9683
if ( cellsizeX < 0 )
9784
{
9885
cellsizeX = -cellsizeX;
9986
}
100-
double cellsizeY = geoTransform[5];
87+
double cellsizeY = mRasterLayer->rasterUnitsPerPixelY();
10188
if ( cellsizeY < 0 )
10289
{
10390
cellsizeY = -cellsizeY;
10491
}
105-
QgsRectangle rasterBBox( geoTransform[0], geoTransform[3] - ( nCellsYGDAL * cellsizeY ),
106-
geoTransform[0] + ( nCellsXGDAL * cellsizeX ), geoTransform[3] );
92+
QgsRectangle rasterBBox = mRasterProvider->extent();
10793

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

226-
227212
//iterate over each polygon
228213
QgsFeatureRequest request;
229214
request.setSubsetOfAttributes( QgsAttributeList() );
@@ -273,22 +258,22 @@ int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
273258
}
274259

275260
//avoid access to cells outside of the raster (may occur because of rounding)
276-
if (( offsetX + nCellsX ) > nCellsXGDAL )
261+
if (( offsetX + nCellsX ) > nCellsXProvider )
277262
{
278-
nCellsX = nCellsXGDAL - offsetX;
263+
nCellsX = nCellsXProvider - offsetX;
279264
}
280-
if (( offsetY + nCellsY ) > nCellsYGDAL )
265+
if (( offsetY + nCellsY ) > nCellsYProvider )
281266
{
282-
nCellsY = nCellsYGDAL - offsetY;
267+
nCellsY = nCellsYProvider - offsetY;
283268
}
284269

285-
statisticsFromMiddlePointTest( rasterBand, featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
270+
statisticsFromMiddlePointTest( featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
286271
rasterBBox, featureStats );
287272

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

@@ -366,7 +351,6 @@ int QgsZonalStatistics::calculateStatistics( QProgressDialog* p )
366351
p->setValue( featureCount );
367352
}
368353

369-
GDALClose( inputDataset );
370354
mPolygonLayer->updateFields();
371355

372356
if ( p && p->wasCanceled() )
@@ -404,12 +388,11 @@ int QgsZonalStatistics::cellInfoForBBox( const QgsRectangle& rasterBBox, const Q
404388
return 0;
405389
}
406390

407-
void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, const QgsGeometry& poly, int pixelOffsetX,
391+
void QgsZonalStatistics::statisticsFromMiddlePointTest( const QgsGeometry& poly, int pixelOffsetX,
408392
int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, FeatureStats &stats )
409393
{
410394
double cellCenterX, cellCenterY;
411395

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

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

416+
QgsRectangle featureBBox = poly.boundingBox().intersect( &rasterBBox );
417+
QgsRectangle intersectBBox = rasterBBox.intersect( &featureBBox );
418+
419+
QgsRasterBlock* block = mRasterProvider->block( mRasterBand, intersectBBox, nCellsX, nCellsY );
433420
for ( int i = 0; i < nCellsY; ++i )
434421
{
435-
if ( GDALRasterIO( band, GF_Read, pixelOffsetX, pixelOffsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 )
436-
!= CPLE_None )
437-
{
438-
continue;
439-
}
440422
cellCenterX = rasterBBox.xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;
441423
for ( int j = 0; j < nCellsX; ++j )
442424
{
443-
if ( validPixel( scanLine[j] ) )
425+
if ( validPixel( block->value( i, j ) ) )
444426
{
445427
GEOSGeom_destroy_r( geosctxt, currentCellCenter );
446428
cellCenterCoords = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
@@ -449,45 +431,46 @@ void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, const QgsGeo
449431
currentCellCenter = GEOSGeom_createPoint_r( geosctxt, cellCenterCoords );
450432
if ( GEOSPreparedContains_r( geosctxt, polyGeosPrepared, currentCellCenter ) )
451433
{
452-
stats.addValue( scanLine[j] );
434+
stats.addValue( block->value( i, j ) );
453435
}
454436
}
455437
cellCenterX += cellSizeX;
456438
}
457439
cellCenterY -= cellSizeY;
458440
}
441+
459442
GEOSGeom_destroy_r( geosctxt, currentCellCenter );
460-
CPLFree( scanLine );
461443
GEOSPreparedGeom_destroy_r( geosctxt, polyGeosPrepared );
462444
GEOSGeom_destroy_r( geosctxt, polyGeos );
445+
delete block;
463446
}
464447

465-
void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, const QgsGeometry& poly, int pixelOffsetX,
448+
void QgsZonalStatistics::statisticsFromPreciseIntersection( const QgsGeometry& poly, int pixelOffsetX,
466449
int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, FeatureStats &stats )
467450
{
468451
stats.reset();
469452

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

474456
double hCellSizeX = cellSizeX / 2.0;
475457
double hCellSizeY = cellSizeY / 2.0;
476458
double pixelArea = cellSizeX * cellSizeY;
477459
double weight = 0;
478460

479-
for ( int row = 0; row < nCellsY; ++row )
461+
QgsRectangle featureBBox = poly.boundingBox().intersect( &rasterBBox );
462+
QgsRectangle intersectBBox = rasterBBox.intersect( &featureBBox );
463+
464+
QgsRasterBlock* block = mRasterProvider->block( mRasterBand, intersectBBox, nCellsX, nCellsY );
465+
for ( int i = 0; i < nCellsY; ++i )
480466
{
481467
double currentX = rasterBBox.xMinimum() + cellSizeX / 2.0 + pixelOffsetX * cellSizeX;
482-
for ( int col = 0; col < nCellsX; ++col )
468+
for ( int j = 0; j < nCellsX; ++j )
483469
{
484-
if ( GDALRasterIO( band, GF_Read, pixelOffsetX + col, pixelOffsetY + row, nCellsX, 1, pixelData, 1, 1, GDT_Float32, 0, 0 ) != CE_None )
470+
if ( !validPixel( block->value( i, j ) ) )
485471
{
486-
QgsDebugMsg( "Raster IO Error" );
487-
}
488-
489-
if ( !validPixel( *pixelData ) )
490472
continue;
473+
}
491474

492475
pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) );
493476
if ( !pixelRectGeometry.isNull() )
@@ -500,7 +483,7 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, const Qg
500483
if ( intersectionArea >= 0.0 )
501484
{
502485
weight = intersectionArea / pixelArea;
503-
stats.addValue( *pixelData, weight );
486+
stats.addValue( block->value( i, j ), weight );
504487
}
505488
}
506489
pixelRectGeometry = QgsGeometry();
@@ -509,7 +492,7 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, const Qg
509492
}
510493
currentY -= cellSizeY;
511494
}
512-
CPLFree( pixelData );
495+
delete block;
513496
}
514497

515498
bool QgsZonalStatistics::validPixel( float value ) const

src/analysis/vector/qgszonalstatistics.h

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,8 @@
2626

2727
class QgsGeometry;
2828
class QgsVectorLayer;
29+
class QgsRasterLayer;
30+
class QgsRasterDataProvider;
2931
class QProgressDialog;
3032
class QgsRectangle;
3133
class QgsField;
@@ -57,7 +59,7 @@ class ANALYSIS_EXPORT QgsZonalStatistics
5759
/**
5860
* Constructor for QgsZonalStatistics.
5961
*/
60-
QgsZonalStatistics( QgsVectorLayer* polygonLayer, const QString& rasterFile, const QString& attributePrefix = "", int rasterBand = 1,
62+
QgsZonalStatistics( QgsVectorLayer* polygonLayer, QgsRasterLayer* rasterLayer, const QString& attributePrefix = "", int rasterBand = 1,
6163
Statistics stats = Statistics( Count | Sum | Mean ) );
6264

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

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

120122
//! Returns statistics with precise pixel - polygon intersection test (slow)
121-
void statisticsFromPreciseIntersection( void* band, const QgsGeometry& poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY,
123+
void statisticsFromPreciseIntersection( const QgsGeometry& poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY,
122124
double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, FeatureStats& stats );
123125

124126
//! Tests whether a pixel's value should be included in the result
125127
bool validPixel( float value ) const;
126128

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

129-
QString mRasterFilePath;
131+
QgsRasterLayer* mRasterLayer;
132+
QgsRasterDataProvider* mRasterProvider;
130133
//! Raster band to calculate statistics from (defaults to 1)
131134
int mRasterBand;
132135
QgsVectorLayer* mPolygonLayer;

tests/src/analysis/testqgszonalstatistics.cpp

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@
1919
#include "qgsapplication.h"
2020
#include "qgsfeatureiterator.h"
2121
#include "qgsvectorlayer.h"
22+
#include "qgsrasterlayer.h"
2223
#include "qgszonalstatistics.h"
2324
#include "qgsproject.h"
2425

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

4344
private:
4445
QgsVectorLayer* mVectorLayer;
45-
QString mRasterPath;
46+
QgsRasterLayer* mRasterLayer;
4647
};
4748

4849
TestQgsZonalStatistics::TestQgsZonalStatistics()
@@ -71,10 +72,9 @@ void TestQgsZonalStatistics::initTestCase()
7172
}
7273

7374
mVectorLayer = new QgsVectorLayer( myTempPath + "polys.shp", QStringLiteral( "poly" ), QStringLiteral( "ogr" ) );
75+
mRasterLayer = new QgsRasterLayer( myTempPath + "edge_problem.asc", QStringLiteral( "raster" ), QStringLiteral( "gdal" ) );
7476
QgsProject::instance()->addMapLayers(
75-
QList<QgsMapLayer *>() << mVectorLayer );
76-
77-
mRasterPath = myTempPath + "edge_problem.asc";
77+
QList<QgsMapLayer *>() << mVectorLayer << mRasterLayer );
7878
}
7979

8080
void TestQgsZonalStatistics::cleanupTestCase()
@@ -84,7 +84,7 @@ void TestQgsZonalStatistics::cleanupTestCase()
8484

8585
void TestQgsZonalStatistics::testStatistics()
8686
{
87-
QgsZonalStatistics zs( mVectorLayer, mRasterPath, QLatin1String( "" ), 1, QgsZonalStatistics::All );
87+
QgsZonalStatistics zs( mVectorLayer, mRasterLayer, QLatin1String( "" ), 1, QgsZonalStatistics::All );
8888
zs.calculateStatistics( nullptr );
8989

9090
QgsFeature f;
@@ -135,7 +135,7 @@ void TestQgsZonalStatistics::testStatistics()
135135
QCOMPARE( f.attribute( "variety" ).toDouble(), 2.0 );
136136

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

141141
request.setFilterFid( 0 );

tests/src/python/test_qgszonalstatistics.py

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@
1515
import qgis # NOQA
1616

1717
from qgis.PyQt.QtCore import QDir, QFile
18-
from qgis.core import QgsVectorLayer, QgsFeature, QgsFeatureRequest
18+
from qgis.core import QgsVectorLayer, QgsRasterLayer, QgsFeature, QgsFeatureRequest
1919
from qgis.analysis import QgsZonalStatistics
2020

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

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

4545
feat = QgsFeature()

0 commit comments

Comments
 (0)