From 278913b7ce526c64125630a76dd21efb204f9f09 Mon Sep 17 00:00:00 2001 From: Alexander Bruy Date: Thu, 26 Jan 2017 16:19:22 +0200 Subject: [PATCH 1/2] use QgsRasterBlock instead of GDAL in zonal statistics (addresses #8736) This should make zonal statistics usable rasters which come from other providers, e.g. WCS. --- python/analysis/vector/qgszonalstatistics.sip | 2 +- .../algs/qgis/ZonalStatisticsQgis.py | 3 +- src/analysis/vector/qgszonalstatistics.cpp | 101 ++++++++---------- src/analysis/vector/qgszonalstatistics.h | 11 +- tests/src/analysis/testqgszonalstatistics.cpp | 12 +-- tests/src/python/test_qgszonalstatistics.py | 6 +- 6 files changed, 61 insertions(+), 74 deletions(-) diff --git a/python/analysis/vector/qgszonalstatistics.sip b/python/analysis/vector/qgszonalstatistics.sip index fbf6d25a910d..9e36ed0c0f58 100644 --- a/python/analysis/vector/qgszonalstatistics.sip +++ b/python/analysis/vector/qgszonalstatistics.sip @@ -30,7 +30,7 @@ class QgsZonalStatistics typedef QFlags 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 diff --git a/python/plugins/processing/algs/qgis/ZonalStatisticsQgis.py b/python/plugins/processing/algs/qgis/ZonalStatisticsQgis.py index ec65dc0c76b6..547ed1e07060 100644 --- a/python/plugins/processing/algs/qgis/ZonalStatisticsQgis.py +++ b/python/plugins/processing/algs/qgis/ZonalStatisticsQgis.py @@ -100,6 +100,7 @@ def processAlgorithm(self, feedback): st = self.getParameterValue(self.STATISTICS) vectorLayer = dataobjects.getObjectFromUri(vectorPath) + rasterLayer = dataobjects.getObjectFromUri(rasterPath) keys = list(self.STATS.keys()) selectedStats = 0 @@ -107,7 +108,7 @@ def processAlgorithm(self, feedback): selectedStats |= self.STATS[keys[i]] zs = QgsZonalStatistics(vectorLayer, - rasterPath, + rasterLayer, columnPrefix, bandNumber, selectedStats) diff --git a/src/analysis/vector/qgszonalstatistics.cpp b/src/analysis/vector/qgszonalstatistics.cpp index 98e7b7dd520d..f80bc0bd5213 100644 --- a/src/analysis/vector/qgszonalstatistics.cpp +++ b/src/analysis/vector/qgszonalstatistics.cpp @@ -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 #include -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 ) @@ -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 ) @@ -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 newFieldList; @@ -223,7 +209,6 @@ int QgsZonalStatistics::calculateStatistics( QProgressDialog* p ) p->setMaximum( featureCount ); } - //iterate over each polygon QgsFeatureRequest request; request.setSubsetOfAttributes( QgsAttributeList() ); @@ -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 ); } @@ -366,7 +351,6 @@ int QgsZonalStatistics::calculateStatistics( QProgressDialog* p ) p->setValue( featureCount ); } - GDALClose( inputDataset ); mPolygonLayer->updateFields(); if ( p && p->wasCanceled() ) @@ -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(); @@ -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 ); @@ -449,26 +431,26 @@ 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; @@ -476,18 +458,19 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, const Qg 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() ) @@ -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(); @@ -509,7 +492,7 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, const Qg } currentY -= cellSizeY; } - CPLFree( pixelData ); + delete block; } bool QgsZonalStatistics::validPixel( float value ) const diff --git a/src/analysis/vector/qgszonalstatistics.h b/src/analysis/vector/qgszonalstatistics.h index 0b0eec60d52f..b623d6dccba3 100644 --- a/src/analysis/vector/qgszonalstatistics.h +++ b/src/analysis/vector/qgszonalstatistics.h @@ -26,6 +26,8 @@ class QgsGeometry; class QgsVectorLayer; +class QgsRasterLayer; +class QgsRasterDataProvider; class QProgressDialog; class QgsRectangle; class QgsField; @@ -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 @@ -114,11 +116,11 @@ 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 @@ -126,7 +128,8 @@ class ANALYSIS_EXPORT QgsZonalStatistics QString getUniqueFieldName( const QString& fieldName, const QList& newFields ); - QString mRasterFilePath; + QgsRasterLayer* mRasterLayer; + QgsRasterDataProvider* mRasterProvider; //! Raster band to calculate statistics from (defaults to 1) int mRasterBand; QgsVectorLayer* mPolygonLayer; diff --git a/tests/src/analysis/testqgszonalstatistics.cpp b/tests/src/analysis/testqgszonalstatistics.cpp index 0f5f7e5d4f6c..0f2889d386bf 100644 --- a/tests/src/analysis/testqgszonalstatistics.cpp +++ b/tests/src/analysis/testqgszonalstatistics.cpp @@ -19,6 +19,7 @@ #include "qgsapplication.h" #include "qgsfeatureiterator.h" #include "qgsvectorlayer.h" +#include "qgsrasterlayer.h" #include "qgszonalstatistics.h" #include "qgsproject.h" @@ -42,7 +43,7 @@ class TestQgsZonalStatistics : public QObject private: QgsVectorLayer* mVectorLayer; - QString mRasterPath; + QgsRasterLayer* mRasterLayer; }; TestQgsZonalStatistics::TestQgsZonalStatistics() @@ -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() << mVectorLayer ); - - mRasterPath = myTempPath + "edge_problem.asc"; + QList() << mVectorLayer << mRasterLayer ); } void TestQgsZonalStatistics::cleanupTestCase() @@ -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; @@ -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 ); diff --git a/tests/src/python/test_qgszonalstatistics.py b/tests/src/python/test_qgszonalstatistics.py index ec22079b35a0..4be3ca32455c 100644 --- a/tests/src/python/test_qgszonalstatistics.py +++ b/tests/src/python/test_qgszonalstatistics.py @@ -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 @@ -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() From d681e9cb060b231db384bd3bafc0c2ce8ba9de13 Mon Sep 17 00:00:00 2001 From: Alexander Bruy Date: Thu, 26 Jan 2017 16:26:53 +0200 Subject: [PATCH 2/2] update docs --- doc/api_break.dox | 3 +++ 1 file changed, 3 insertions(+) diff --git a/doc/api_break.dox b/doc/api_break.dox index 0d50e909e8c7..f435e708d7ec 100644 --- a/doc/api_break.dox +++ b/doc/api_break.dox @@ -2014,6 +2014,9 @@ Triangulation {#qgis_api_break_3_0_Triangulation} - 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. +QgsZonalStatistics {#qgis_api_break_3_0_QgsZonalStatistics} +------------------ +- QgsZonalStatistics() сonstructor now accepts pointer to the QgsRasterLayer instance instead of path to the raster file QGIS 2.6 {#qgis_api_break_2_6}