QGIS API Documentation 3.41.0-Master (af5edcb665c)
Loading...
Searching...
No Matches
qgsalgorithmrasterzonalstats.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsalgorithmrasterzonalstats.cpp
3 ---------------------
4 begin : December 2018
5 copyright : (C) 2018 by Nyall Dawson
6 email : nyall dot dawson at gmail dot com
7 ***************************************************************************/
8
9/***************************************************************************
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 ***************************************************************************/
17
19#include "qgsstringutils.h"
21#include "qgsrasterprojector.h"
22#include "qgsunittypes.h"
23
25
26QString QgsRasterLayerZonalStatsAlgorithm::name() const
27{
28 return QStringLiteral( "rasterlayerzonalstats" );
29}
30
31QString QgsRasterLayerZonalStatsAlgorithm::displayName() const
32{
33 return QObject::tr( "Raster layer zonal statistics" );
34}
35
36QStringList QgsRasterLayerZonalStatsAlgorithm::tags() const
37{
38 return QObject::tr( "count,area,statistics,stats,zones,categories,minimum,maximum,mean,sum,total" ).split( ',' );
39}
40
41QString QgsRasterLayerZonalStatsAlgorithm::group() const
42{
43 return QObject::tr( "Raster analysis" );
44}
45
46QString QgsRasterLayerZonalStatsAlgorithm::groupId() const
47{
48 return QStringLiteral( "rasteranalysis" );
49}
50
51void QgsRasterLayerZonalStatsAlgorithm::initAlgorithm( const QVariantMap & )
52{
53 addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT" ), QObject::tr( "Input layer" ) ) );
54 addParameter( new QgsProcessingParameterBand( QStringLiteral( "BAND" ), QObject::tr( "Band number" ), 1, QStringLiteral( "INPUT" ) ) );
55 addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "ZONES" ), QObject::tr( "Zones layer" ) ) );
56 addParameter( new QgsProcessingParameterBand( QStringLiteral( "ZONES_BAND" ), QObject::tr( "Zones band number" ), 1, QStringLiteral( "ZONES" ) ) );
57
58 std::unique_ptr<QgsProcessingParameterEnum> refParam = std::make_unique<QgsProcessingParameterEnum>( QStringLiteral( "REF_LAYER" ), QObject::tr( "Reference layer" ), QStringList() << QObject::tr( "Input layer" ) << QObject::tr( "Zones layer" ), false, 0 );
59 refParam->setFlags( refParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
60 addParameter( refParam.release() );
61
62 addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT_TABLE" ), QObject::tr( "Statistics" ), Qgis::ProcessingSourceType::Vector ) );
63
64 addOutput( new QgsProcessingOutputString( QStringLiteral( "EXTENT" ), QObject::tr( "Extent" ) ) );
65 addOutput( new QgsProcessingOutputString( QStringLiteral( "CRS_AUTHID" ), QObject::tr( "CRS authority identifier" ) ) );
66 addOutput( new QgsProcessingOutputNumber( QStringLiteral( "WIDTH_IN_PIXELS" ), QObject::tr( "Width in pixels" ) ) );
67 addOutput( new QgsProcessingOutputNumber( QStringLiteral( "HEIGHT_IN_PIXELS" ), QObject::tr( "Height in pixels" ) ) );
68 addOutput( new QgsProcessingOutputNumber( QStringLiteral( "TOTAL_PIXEL_COUNT" ), QObject::tr( "Total pixel count" ) ) );
69 addOutput( new QgsProcessingOutputNumber( QStringLiteral( "NODATA_PIXEL_COUNT" ), QObject::tr( "NoData pixel count" ) ) );
70}
71
72QString QgsRasterLayerZonalStatsAlgorithm::shortDescription() const
73{
74 return QObject::tr( "Calculates statistics for a raster layer's values, categorized by zones defined in another raster layer." );
75}
76
77QString QgsRasterLayerZonalStatsAlgorithm::shortHelpString() const
78{
79 return QObject::tr( "This algorithm calculates statistics for a raster layer's values, categorized by zones defined in another raster layer.\n\n"
80 "If the reference layer parameter is set to \"Input layer\", then zones are determined by sampling the zone raster layer value at the centroid of each pixel from the source raster layer.\n\n"
81 "If the reference layer parameter is set to \"Zones layer\", then the input raster layer will be sampled at the centroid of each pixel from the zones raster layer.\n\n"
82 "If either the source raster layer or the zone raster layer value is NoData for a pixel, that pixel's value will be skipped and not included in the calculated statistics." );
83}
84
85QgsRasterLayerZonalStatsAlgorithm *QgsRasterLayerZonalStatsAlgorithm::createInstance() const
86{
87 return new QgsRasterLayerZonalStatsAlgorithm();
88}
89
90bool QgsRasterLayerZonalStatsAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
91{
92 mRefLayer = static_cast<RefLayer>( parameterAsEnum( parameters, QStringLiteral( "REF_LAYER" ), context ) );
93
94 QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT" ), context );
95 const int band = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
96
97 if ( !layer )
98 throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT" ) ) );
99
100 mBand = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
101 if ( mBand < 1 || mBand > layer->bandCount() )
102 throw QgsProcessingException( QObject::tr( "Invalid band number for BAND (%1): Valid values for input raster are 1 to %2" ).arg( mBand ).arg( layer->bandCount() ) );
103
104 mHasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( band );
105
106 QgsRasterLayer *zonesLayer = parameterAsRasterLayer( parameters, QStringLiteral( "ZONES" ), context );
107
108 if ( !zonesLayer )
109 throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "ZONES" ) ) );
110
111 mZonesBand = parameterAsInt( parameters, QStringLiteral( "ZONES_BAND" ), context );
112 if ( mZonesBand < 1 || mZonesBand > zonesLayer->bandCount() )
113 throw QgsProcessingException( QObject::tr( "Invalid band number for ZONES_BAND (%1): Valid values for input raster are 1 to %2" ).arg( mZonesBand ).arg( zonesLayer->bandCount() ) );
114 mZonesHasNoDataValue = zonesLayer->dataProvider()->sourceHasNoDataValue( band );
115
116 mSourceDataProvider.reset( layer->dataProvider()->clone() );
117 mSourceInterface = mSourceDataProvider.get();
118 mZonesDataProvider.reset( zonesLayer->dataProvider()->clone() );
119 mZonesInterface = mZonesDataProvider.get();
120
121 switch ( mRefLayer )
122 {
123 case Source:
124 mCrs = layer->crs();
125 mRasterUnitsPerPixelX = layer->rasterUnitsPerPixelX();
126 mRasterUnitsPerPixelY = layer->rasterUnitsPerPixelY();
127 mLayerWidth = layer->width();
128 mLayerHeight = layer->height();
129 mExtent = layer->extent();
130
131 // add projector if necessary
132 if ( layer->crs() != zonesLayer->crs() )
133 {
134 mProjector = std::make_unique<QgsRasterProjector>();
135 mProjector->setInput( mZonesDataProvider.get() );
136 mProjector->setCrs( zonesLayer->crs(), layer->crs(), context.transformContext() );
137 mZonesInterface = mProjector.get();
138 }
139 break;
140
141 case Zones:
142 mCrs = zonesLayer->crs();
143 mRasterUnitsPerPixelX = zonesLayer->rasterUnitsPerPixelX();
144 mRasterUnitsPerPixelY = zonesLayer->rasterUnitsPerPixelY();
145 mLayerWidth = zonesLayer->width();
146 mLayerHeight = zonesLayer->height();
147 mExtent = zonesLayer->extent();
148
149 // add projector if necessary
150 if ( layer->crs() != zonesLayer->crs() )
151 {
152 mProjector = std::make_unique<QgsRasterProjector>();
153 mProjector->setInput( mSourceDataProvider.get() );
154 mProjector->setCrs( layer->crs(), zonesLayer->crs(), context.transformContext() );
155 mSourceInterface = mProjector.get();
156 }
157 break;
158 }
159
160 return true;
161}
162
163QVariantMap QgsRasterLayerZonalStatsAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
164{
165 QString areaUnit = QgsUnitTypes::toAbbreviatedString( QgsUnitTypes::distanceToAreaUnit( mCrs.mapUnits() ) );
166
167 QString tableDest;
168 std::unique_ptr<QgsFeatureSink> sink;
169 if ( parameters.contains( QStringLiteral( "OUTPUT_TABLE" ) ) && parameters.value( QStringLiteral( "OUTPUT_TABLE" ) ).isValid() )
170 {
171 QgsFields outFields;
172 outFields.append( QgsField( QStringLiteral( "zone" ), QMetaType::Type::Double, QString(), 20, 8 ) );
173 outFields.append( QgsField( areaUnit.isEmpty() ? "area" : areaUnit.replace( QStringLiteral( "²" ), QStringLiteral( "2" ) ), QMetaType::Type::Double, QString(), 20, 8 ) );
174 outFields.append( QgsField( QStringLiteral( "sum" ), QMetaType::Type::Double, QString(), 20, 8 ) );
175 outFields.append( QgsField( QStringLiteral( "count" ), QMetaType::Type::LongLong, QString(), 20 ) );
176 outFields.append( QgsField( QStringLiteral( "min" ), QMetaType::Type::Double, QString(), 20, 8 ) );
177 outFields.append( QgsField( QStringLiteral( "max" ), QMetaType::Type::Double, QString(), 20, 8 ) );
178 outFields.append( QgsField( QStringLiteral( "mean" ), QMetaType::Type::Double, QString(), 20, 8 ) );
179
180 sink.reset( parameterAsSink( parameters, QStringLiteral( "OUTPUT_TABLE" ), context, tableDest, outFields, Qgis::WkbType::NoGeometry, QgsCoordinateReferenceSystem() ) );
181 if ( !sink )
182 throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT_TABLE" ) ) );
183 }
184
185 struct StatCalculator
186 {
187 // only calculate cheap stats-- we cannot calculate stats which require holding values in memory -- because otherwise we'll end
188 // up trying to store EVERY pixel value from the input in memory
190 };
191 QHash<double, StatCalculator> zoneStats;
192 qgssize noDataCount = 0;
193
194 const qgssize layerSize = static_cast<qgssize>( mLayerWidth ) * static_cast<qgssize>( mLayerHeight );
197 const int nbBlocksWidth = static_cast<int>( std::ceil( 1.0 * mLayerWidth / maxWidth ) );
198 const int nbBlocksHeight = static_cast<int>( std::ceil( 1.0 * mLayerHeight / maxHeight ) );
199 const int nbBlocks = nbBlocksWidth * nbBlocksHeight;
200
201 QgsRasterIterator iter = mRefLayer == Source ? QgsRasterIterator( mSourceInterface )
202 : QgsRasterIterator( mZonesInterface );
203 iter.startRasterRead( mRefLayer == Source ? mBand : mZonesBand, mLayerWidth, mLayerHeight, mExtent );
204
205 int iterLeft = 0;
206 int iterTop = 0;
207 int iterCols = 0;
208 int iterRows = 0;
209 QgsRectangle blockExtent;
210 std::unique_ptr<QgsRasterBlock> rasterBlock;
211 std::unique_ptr<QgsRasterBlock> zonesRasterBlock;
212 bool isNoData = false;
213 while ( true )
214 {
215 if ( mRefLayer == Source )
216 {
217 if ( !iter.readNextRasterPart( mBand, iterCols, iterRows, rasterBlock, iterLeft, iterTop, &blockExtent ) )
218 break;
219
220 zonesRasterBlock.reset( mZonesInterface->block( mZonesBand, blockExtent, iterCols, iterRows ) );
221 }
222 else
223 {
224 if ( !iter.readNextRasterPart( mZonesBand, iterCols, iterRows, zonesRasterBlock, iterLeft, iterTop, &blockExtent ) )
225 break;
226
227 rasterBlock.reset( mSourceInterface->block( mBand, blockExtent, iterCols, iterRows ) );
228 }
229 if ( !zonesRasterBlock || !rasterBlock )
230 continue;
231
232 feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
233 if ( !rasterBlock->isValid() || rasterBlock->isEmpty() || !zonesRasterBlock->isValid() || zonesRasterBlock->isEmpty() )
234 continue;
235
236 for ( int row = 0; row < iterRows; row++ )
237 {
238 if ( feedback->isCanceled() )
239 break;
240
241 for ( int column = 0; column < iterCols; column++ )
242 {
243 const double value = rasterBlock->valueAndNoData( row, column, isNoData );
244 if ( mHasNoDataValue && isNoData )
245 {
246 noDataCount += 1;
247 continue;
248 }
249 const double zone = zonesRasterBlock->valueAndNoData( row, column, isNoData );
250 if ( mZonesHasNoDataValue && isNoData )
251 {
252 noDataCount += 1;
253 continue;
254 }
255 zoneStats[zone].s.addValue( value );
256 }
257 }
258 }
259
260 QVariantMap outputs;
261 outputs.insert( QStringLiteral( "EXTENT" ), mExtent.toString() );
262 outputs.insert( QStringLiteral( "CRS_AUTHID" ), mCrs.authid() );
263 outputs.insert( QStringLiteral( "WIDTH_IN_PIXELS" ), mLayerWidth );
264 outputs.insert( QStringLiteral( "HEIGHT_IN_PIXELS" ), mLayerHeight );
265 outputs.insert( QStringLiteral( "TOTAL_PIXEL_COUNT" ), layerSize );
266 outputs.insert( QStringLiteral( "NODATA_PIXEL_COUNT" ), noDataCount );
267
268 const double pixelArea = mRasterUnitsPerPixelX * mRasterUnitsPerPixelY;
269
270 for ( auto it = zoneStats.begin(); it != zoneStats.end(); ++it )
271 {
272 QgsFeature f;
273 it->s.finalize();
274 f.setAttributes( QgsAttributes() << it.key() << it->s.count() * pixelArea << it->s.sum() << it->s.count() << it->s.min() << it->s.max() << it->s.mean() );
275 if ( !sink->addFeature( f, QgsFeatureSink::FastInsert ) )
276 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT_TABLE" ) ) );
277 sink->finalize();
278 }
279 outputs.insert( QStringLiteral( "OUTPUT_TABLE" ), tableDest );
280
281 return outputs;
282}
283
284
@ Vector
Tables (i.e. vector layers with or without geometry). When used for a sink this indicates the sink ha...
@ Mean
Mean of values.
@ Max
Max of values.
@ Min
Min of values.
@ Sum
Sum of values.
@ NoGeometry
No geometry.
@ Advanced
Parameter is an advanced parameter which should be hidden from users by default.
A vector of attributes.
This class represents a coordinate reference system (CRS).
@ FastInsert
Use faster inserts, at the cost of updating the passed features to reflect changes made at the provid...
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition qgsfeature.h:58
void setAttributes(const QgsAttributes &attrs)
Sets the feature's attributes.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition qgsfeedback.h:53
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition qgsfeedback.h:61
Encapsulate a field in an attribute table or data source.
Definition qgsfield.h:53
Container of fields for a vector layer.
Definition qgsfields.h:46
bool append(const QgsField &field, Qgis::FieldOrigin origin=Qgis::FieldOrigin::Provider, int originIndex=-1)
Appends a field.
Definition qgsfields.cpp:70
virtual QgsRectangle extent() const
Returns the extent of the layer.
QgsCoordinateReferenceSystem crs
Definition qgsmaplayer.h:83
Contains information about the context in which a processing algorithm is executed.
QgsCoordinateTransformContext transformContext() const
Returns the coordinate transform context.
Custom exception class for processing related exceptions.
Base class for providing feedback from a processing algorithm.
A numeric output for processing algorithms.
A string output for processing algorithms.
A raster band parameter for Processing algorithms.
A feature sink output for processing algorithms.
A raster layer parameter for processing algorithms.
QgsRasterDataProvider * clone() const override=0
Clone itself, create deep copy.
virtual bool sourceHasNoDataValue(int bandNo) const
Returns true if source band has no data value.
Iterator for sequentially processing raster cells.
static const int DEFAULT_MAXIMUM_TILE_WIDTH
Default maximum tile width.
bool readNextRasterPart(int bandNumber, int &nCols, int &nRows, QgsRasterBlock **block, int &topLeftCol, int &topLeftRow)
Fetches next part of raster data, caller takes ownership of the block and caller should delete the bl...
static const int DEFAULT_MAXIMUM_TILE_HEIGHT
Default maximum tile height.
void startRasterRead(int bandNumber, qgssize nCols, qgssize nRows, const QgsRectangle &extent, QgsRasterBlockFeedback *feedback=nullptr)
Start reading of raster band.
Represents a raster layer.
int height() const
Returns the height of the (unclipped) raster.
int bandCount() const
Returns the number of bands in this layer.
double rasterUnitsPerPixelX() const
Returns the number of raster units per each raster pixel in X axis.
QgsRasterDataProvider * dataProvider() override
Returns the source data provider.
double rasterUnitsPerPixelY() const
Returns the number of raster units per each raster pixel in Y axis.
int width() const
Returns the width of the (unclipped) raster.
A rectangle specified with double values.
Calculator for summary statistics for a list of doubles.
static Q_INVOKABLE QString toAbbreviatedString(Qgis::DistanceUnit unit)
Returns a translated abbreviation representing a distance unit.
static Q_INVOKABLE Qgis::AreaUnit distanceToAreaUnit(Qgis::DistanceUnit distanceUnit)
Converts a distance unit to its corresponding area unit, e.g., meters to square meters.
unsigned long long qgssize
Qgssize is used instead of size_t, because size_t is stdlib type, unknown by SIP, and it would be har...
Definition qgis.h:6614