Add thread safe method of constructing QgsZonalStatistics,

using a QgsRasterInterface instead of a QgsRasterLayer

python/analysis/auto_generated/vector/qgszonalstatistics.sip.in

@@ -48,9 +48,36 @@ A class that calculates raster statistics (count, sum, mean) for a polygon or mu
                         int rasterBand = 1,
                         QgsZonalStatistics::Statistics stats = QgsZonalStatistics::Statistics( QgsZonalStatistics::Count | QgsZonalStatistics::Sum | QgsZonalStatistics::Mean ) );
 %Docstring
-Constructor for QgsZonalStatistics.
+Convenience constructor for QgsZonalStatistics, using an input raster layer.
+
+The raster layer must exist for the lifetime of the zonal statistics calculation.
+
+.. warning::
+
+   Constructing QgsZonalStatistics using this method is not thread safe, and
+   the constructor which accepts a QgsRasterInterface should be used instead.
 %End
 
+    QgsZonalStatistics( QgsVectorLayer *polygonLayer,
+                        QgsRasterInterface *rasterInterface,
+                        const QgsCoordinateReferenceSystem &rasterCrs,
+                        double rasterUnitsPerPixelX,
+                        double rasterUnitsPerPixelY,
+                        const QString &attributePrefix = QString(),
+                        int rasterBand = 1,
+                        QgsZonalStatistics::Statistics stats = QgsZonalStatistics::Statistics( QgsZonalStatistics::Count | QgsZonalStatistics::Sum | QgsZonalStatistics::Mean ) );
+%Docstring
+Constructor for QgsZonalStatistics, using a QgsRasterInterface.
+
+.. warning::
+
+   The raster interface must exist for the lifetime of the zonal statistics calculation. For thread
+   safe use, always use a cloned raster interface.
+
+.. versionadded:: 3.2
+%End
+
+
     int calculateStatistics( QgsFeedback *feedback );
 %Docstring
 Starts the calculation

src/analysis/vector/qgszonalstatistics.cpp

@@ -31,12 +31,37 @@
 #include <QFile>
 
 QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer *polygonLayer, QgsRasterLayer *rasterLayer, const QString &attributePrefix, int rasterBand, QgsZonalStatistics::Statistics stats )
-  : mRasterLayer( rasterLayer )
+  : QgsZonalStatistics( polygonLayer,
+                        rasterLayer ? rasterLayer->dataProvider() : nullptr,
+                        rasterLayer ? rasterLayer->crs() : QgsCoordinateReferenceSystem(),
+                        rasterLayer ? rasterLayer->rasterUnitsPerPixelX() : 0,
+                        rasterLayer ? rasterLayer->rasterUnitsPerPixelY() : 0,
+                        attributePrefix,
+                        rasterBand,
+                        stats )
+{
+}
+
+QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer *polygonLayer, QgsRasterInterface *rasterInterface,
+                                        const QgsCoordinateReferenceSystem &rasterCrs, double rasterUnitsPerPixelX, double rasterUnitsPerPixelY, const QString &attributePrefix, int rasterBand, QgsZonalStatistics::Statistics stats )
+  : mRasterInterface( rasterInterface )
+  , mRasterCrs( rasterCrs )
+  , mCellSizeX( rasterUnitsPerPixelX )
+  , mCellSizeY( rasterUnitsPerPixelY )
   , mRasterBand( rasterBand )
   , mPolygonLayer( polygonLayer )
   , mAttributePrefix( attributePrefix )
   , mStatistics( stats )
-{}
+{
+  if ( mCellSizeX < 0 )
+  {
+    mCellSizeX = -mCellSizeX;
+  }
+  if ( mCellSizeY < 0 )
+  {
+    mCellSizeY = -mCellSizeY;
+  }
+}
 
 int QgsZonalStatistics::calculateStatistics( QgsFeedback *feedback )
 {
@@ -51,32 +76,21 @@ int QgsZonalStatistics::calculateStatistics( QgsFeedback *feedback )
     return 2;
   }
 
-  if ( !mRasterLayer )
+  if ( !mRasterInterface )
   {
     return 3;
   }
 
-  if ( mRasterLayer->bandCount() < mRasterBand )
+  if ( mRasterInterface->bandCount() < mRasterBand )
   {
     return 4;
   }
 
-  mRasterProvider = mRasterLayer->dataProvider();
-
   //get geometry info about raster layer
-  int nCellsXProvider = mRasterProvider->xSize();
+  int nCellsXProvider = mRasterInterface->xSize();
-  int nCellsYProvider = mRasterProvider->ySize();
+  int nCellsYProvider = mRasterInterface->ySize();
-  double cellsizeX = mRasterLayer->rasterUnitsPerPixelX();
+
-  if ( cellsizeX < 0 )
+  QgsRectangle rasterBBox = mRasterInterface->extent();
-  {
-    cellsizeX = -cellsizeX;
-  }
-  double cellsizeY = mRasterLayer->rasterUnitsPerPixelY();
-  if ( cellsizeY < 0 )
-  {
-    cellsizeY = -cellsizeY;
-  }
-  QgsRectangle rasterBBox = mRasterProvider->extent();
 
   //add the new fields to the provider
   QList<QgsField> newFieldList;
@@ -204,7 +218,7 @@ int QgsZonalStatistics::calculateStatistics( QgsFeedback *feedback )
   //iterate over each polygon
   QgsFeatureRequest request;
   request.setSubsetOfAttributes( QgsAttributeList() );
-  request.setDestinationCrs( mRasterLayer->crs(), QgsProject::instance()->transformContext() );
+  request.setDestinationCrs( mRasterCrs, QgsProject::instance()->transformContext() );
   QgsFeatureIterator fi = vectorProvider->getFeatures( request );
   QgsFeature f;
 
@@ -245,15 +259,15 @@ int QgsZonalStatistics::calculateStatistics( QgsFeedback *feedback )
     }
 
     int offsetX, offsetY, nCellsX, nCellsY;
-    cellInfoForBBox( rasterBBox, featureRect, cellsizeX, cellsizeY, offsetX, offsetY, nCellsX, nCellsY, nCellsXProvider, nCellsYProvider );
+    cellInfoForBBox( rasterBBox, featureRect, mCellSizeX, mCellSizeY, offsetX, offsetY, nCellsX, nCellsY, nCellsXProvider, nCellsYProvider );
 
-    statisticsFromMiddlePointTest( featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
+    statisticsFromMiddlePointTest( featureGeometry, offsetX, offsetY, nCellsX, nCellsY, mCellSizeX, mCellSizeY,
                                    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( featureGeometry, offsetX, offsetY, nCellsX, nCellsY, cellsizeX, cellsizeY,
+      statisticsFromPreciseIntersection( featureGeometry, offsetX, offsetY, nCellsX, nCellsY, mCellSizeX, mCellSizeY,
                                          rasterBBox, featureStats );
     }
 
@@ -394,7 +408,7 @@ void QgsZonalStatistics::statisticsFromMiddlePointTest( const QgsGeometry &poly,
   QgsRectangle featureBBox = poly.boundingBox().intersect( &rasterBBox );
   QgsRectangle intersectBBox = rasterBBox.intersect( &featureBBox );
 
-  std::unique_ptr< QgsRasterBlock > block( mRasterProvider->block( mRasterBand, intersectBBox, nCellsX, nCellsY ) );
+  std::unique_ptr< QgsRasterBlock > block( mRasterInterface->block( mRasterBand, intersectBBox, nCellsX, nCellsY ) );
   for ( int i = 0; i < nCellsY; ++i )
   {
     cellCenterX = rasterBBox.xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;
@@ -438,7 +452,7 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( const QgsGeometry &p
   }
   polyEngine->prepareGeometry();
 
-  std::unique_ptr< QgsRasterBlock > block( mRasterProvider->block( mRasterBand, intersectBBox, nCellsX, nCellsY ) );
+  std::unique_ptr< QgsRasterBlock > block( mRasterInterface->block( mRasterBand, intersectBBox, nCellsX, nCellsY ) );
   for ( int i = 0; i < nCellsY; ++i )
   {
     double currentX = rasterBBox.xMinimum() + cellSizeX / 2.0 + pixelOffsetX * cellSizeX;

src/analysis/vector/qgszonalstatistics.h

@@ -26,10 +26,12 @@
 
 #include "qgis_analysis.h"
 #include "qgsfeedback.h"
+#include "qgscoordinatereferencesystem.h"
 
 class QgsGeometry;
 class QgsVectorLayer;
 class QgsRasterLayer;
+class QgsRasterInterface;
 class QgsRasterDataProvider;
 class QgsRectangle;
 class QgsField;
@@ -61,7 +63,12 @@ class ANALYSIS_EXPORT QgsZonalStatistics
     Q_DECLARE_FLAGS( Statistics, Statistic )
 
     /**
-     * Constructor for QgsZonalStatistics.
+     * Convenience constructor for QgsZonalStatistics, using an input raster layer.
+     *
+     * The raster layer must exist for the lifetime of the zonal statistics calculation.
+     *
+     * \warning Constructing QgsZonalStatistics using this method is not thread safe, and
+     * the constructor which accepts a QgsRasterInterface should be used instead.
      */
     QgsZonalStatistics( QgsVectorLayer *polygonLayer,
                         QgsRasterLayer *rasterLayer,
@@ -69,6 +76,24 @@ class ANALYSIS_EXPORT QgsZonalStatistics
                         int rasterBand = 1,
                         QgsZonalStatistics::Statistics stats = QgsZonalStatistics::Statistics( QgsZonalStatistics::Count | QgsZonalStatistics::Sum | QgsZonalStatistics::Mean ) );
 
+    /**
+     * Constructor for QgsZonalStatistics, using a QgsRasterInterface.
+     *
+     * \warning The raster interface must exist for the lifetime of the zonal statistics calculation. For thread
+     * safe use, always use a cloned raster interface.
+     *
+     * \since QGIS 3.2
+     */
+    QgsZonalStatistics( QgsVectorLayer *polygonLayer,
+                        QgsRasterInterface *rasterInterface,
+                        const QgsCoordinateReferenceSystem &rasterCrs,
+                        double rasterUnitsPerPixelX,
+                        double rasterUnitsPerPixelY,
+                        const QString &attributePrefix = QString(),
+                        int rasterBand = 1,
+                        QgsZonalStatistics::Statistics stats = QgsZonalStatistics::Statistics( QgsZonalStatistics::Count | QgsZonalStatistics::Sum | QgsZonalStatistics::Mean ) );
+
+
     /**
      * Starts the calculation
       \returns 0 in case of success*/
@@ -147,8 +172,12 @@ class ANALYSIS_EXPORT QgsZonalStatistics
 
     QString getUniqueFieldName( const QString &fieldName, const QList<QgsField> &newFields );
 
-    QgsRasterLayer *mRasterLayer = nullptr;
+    QgsRasterInterface *mRasterInterface = nullptr;
-    QgsRasterDataProvider *mRasterProvider = nullptr;
+    QgsCoordinateReferenceSystem mRasterCrs;
+
+    double mCellSizeX = 0;
+    double mCellSizeY = 0;
+
     //! Raster band to calculate statistics
     int mRasterBand = 0;
     QgsVectorLayer *mPolygonLayer = nullptr;