Skip to content
Permalink
Browse files
Fix zonal statistics does not correctly handle coordinate transforms
Fixes #26858
  • Loading branch information
nyalldawson committed Jun 13, 2021
1 parent eb1a55a commit c82678dd6136de79cd76206c9ca25b1f9e34abd2
@@ -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 ] ] ] } }
]
}
@@ -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 ] ] ] } }
]
}
@@ -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:
@@ -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() );
@@ -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;

0 comments on commit c82678d

Please sign in to comment.