Use native QGIS raster API instead of GDAL API in zonal statistics

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. <!--#spellok-->
 
+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}

python/analysis/vector/qgszonalstatistics.sip

@@ -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

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

@@ -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)

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 <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 )
@@ -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<QgsField> 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,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() )
@@ -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

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,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;

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<QgsMapLayer *>() << mVectorLayer );
-
-  mRasterPath = myTempPath + "edge_problem.asc";
+    QList<QgsMapLayer *>() << 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 );

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()