Fix zonal statistics does not correctly handle coordinate transforms

Fixes #26858
qgis · Jun 10, 2021 · 3e5f9ef · 3e5f9ef
1 parent 607b44b
commit 3e5f9ef
Show file tree

Hide file tree

Showing 5 changed files with 67 additions and 2 deletions.
diff --git a/python/plugins/processing/tests/testdata/custom/dem_zones.geojson b/python/plugins/processing/tests/testdata/custom/dem_zones.geojson
@@ -0,0 +1,9 @@
+{
+"type": "FeatureCollection",
+"name": "dem_zones",
+"features": [
+{ "type": "Feature", "properties": { "zone": "zone 1" }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 18.670940771920002, 45.789720648616203 ], [ 18.670601888719101, 45.7786128103632 ], [ 18.680241233101299, 45.778386888229299 ], [ 18.680580116302298, 45.790134839195098 ], [ 18.670940771920002, 45.789720648616203 ] ] ] } },
+{ "type": "Feature", "properties": { "zone": "zone 2" }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 18.696733548880299, 45.804631509457501 ], [ 18.6881861614789, 45.8046691631464 ], [ 18.6969218173253, 45.787687349410596 ], [ 18.696733548880299, 45.804631509457501 ] ] ] } },
+{ "type": "Feature", "properties": { "zone": "zone 3" }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 18.688731551630401, 45.791979281616896 ], [ 18.688842747680699, 45.791921624405603 ], [ 18.6888309809029, 45.791878087327703 ], [ 18.6887509668137, 45.791858672144301 ], [ 18.688731551630401, 45.791979281616896 ] ] ] } }
+]
+}
diff --git a/python/plugins/processing/tests/testdata/expected/zonal_stats_crs.geojson b/python/plugins/processing/tests/testdata/expected/zonal_stats_crs.geojson
@@ -0,0 +1,10 @@
+{
+"type": "FeatureCollection",
+"name": "zonal_stats_crs",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
+"features": [
+{ "type": "Feature", "properties": { "zone": "zone 1", "_count": 11016.0, "_sum": 1717611.2636108398, "_mean": 155.91968623918299 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 18.670940771920002, 45.789720648616203 ], [ 18.670601888719101, 45.7786128103632 ], [ 18.680241233101299, 45.778386888229299 ], [ 18.680580116302298, 45.790134839195098 ], [ 18.670940771920002, 45.789720648616203 ] ] ] } },
+{ "type": "Feature", "properties": { "zone": "zone 2", "_count": 7239.0, "_sum": 1348285.545715332, "_mean": 186.25301087378534 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 18.696733548880299, 45.804631509457501 ], [ 18.6881861614789, 45.8046691631464 ], [ 18.6969218173253, 45.787687349410596 ], [ 18.696733548880299, 45.804631509457501 ] ] ] } },
+{ "type": "Feature", "properties": { "zone": "zone 3", "_count": 0.77731787596712487, "_sum": 128.34554148638338, "_mean": 165.11332809206039 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 18.688731551630401, 45.791979281616896 ], [ 18.688842747680699, 45.791921624405603 ], [ 18.6888309809029, 45.791878087327703 ], [ 18.6887509668137, 45.791858672144301 ], [ 18.688731551630401, 45.791979281616896 ] ] ] } }
+]
+}
diff --git a/python/plugins/processing/tests/testdata/qgis_algorithm_tests3.yaml b/python/plugins/processing/tests/testdata/qgis_algorithm_tests3.yaml
@@ -118,6 +118,34 @@ tests:
           geometry:
             precision: 5
 
+  - algorithm: native:zonalstatisticsfb
+    name: Zonal stats (feature based, with transform)
+    params:
+      COLUMN_PREFIX: _
+      INPUT:
+        name: custom/dem_zones.geojson
+        type: vector
+      INPUT_RASTER:
+        name: custom/dem_crs.tif
+        type: raster
+      RASTER_BAND: 1
+      STATISTICS:
+      - 0
+      - 1
+      - 2
+    results:
+      OUTPUT:
+        name: expected/zonal_stats_crs.geojson
+        type: vector
+        compare:
+          fields:
+            _count:
+              precision: 3
+            _sum:
+              precision: 3
+            _mean:
+              precision: 3
+
   - algorithm: qgis:zonalhistogram
     name: Zonal histogram
     params:

diff --git a/src/analysis/processing/qgsalgorithmzonalstatisticsfeaturebased.cpp b/src/analysis/processing/qgsalgorithmzonalstatisticsfeaturebased.cpp
@@ -158,12 +158,28 @@ bool QgsZonalStatisticsFeatureBasedAlgorithm::prepareAlgorithm( const QVariantMa
 
 QgsFeatureList QgsZonalStatisticsFeatureBasedAlgorithm::processFeature( const QgsFeature &feature, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
 {
-  Q_UNUSED( context )
+  if ( !mCreatedTransform )
+  {
+    mCreatedTransform = true;
+    mFeatureToRasterTransform = QgsCoordinateTransform( sourceCrs(), mCrs, context.transformContext() );
+  }
+
   Q_UNUSED( feedback )
   QgsAttributes attributes = feature.attributes();
   attributes.resize( mOutputFields.size() );
 
-  QMap<QgsZonalStatistics::Statistic, QVariant> results = QgsZonalStatistics::calculateStatistics( mRaster.get(), feature.geometry(), mPixelSizeX, mPixelSizeY, mBand, mStats );
+  QgsGeometry geometry = feature.geometry();
+  try
+  {
+    geometry.transform( mFeatureToRasterTransform );
+  }
+  catch ( QgsCsException & )
+  {
+    if ( feedback )
+      feedback->reportError( QObject::tr( "Encountered a transform error when reprojecting feature with id %1." ).arg( feature.id() ) );
+  }
+
+  QMap<QgsZonalStatistics::Statistic, QVariant> results = QgsZonalStatistics::calculateStatistics( mRaster.get(), geometry, mPixelSizeX, mPixelSizeY, mBand, mStats );
   for ( auto result = results.constBegin(); result != results.constEnd(); ++result )
   {
     attributes.replace( mStatFieldsMapping.value( result.key() ), result.value() );

diff --git a/src/analysis/processing/qgsalgorithmzonalstatisticsfeaturebased.h b/src/analysis/processing/qgsalgorithmzonalstatisticsfeaturebased.h
@@ -62,6 +62,8 @@ class QgsZonalStatisticsFeatureBasedAlgorithm : public QgsProcessingFeatureBased
     QString mPrefix;
     QgsZonalStatistics::Statistics mStats = QgsZonalStatistics::All;
     QgsCoordinateReferenceSystem mCrs;
+    bool mCreatedTransform = false;
+    QgsCoordinateTransform mFeatureToRasterTransform;
     double mPixelSizeX;
     double mPixelSizeY;
     QgsFields mOutputFields;