QGIS API Documentation 3.39.0-Master (67e056379ed)
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" ),
54 QObject::tr( "Input layer" ) ) );
55 addParameter( new QgsProcessingParameterBand( QStringLiteral( "BAND" ),
56 QObject::tr( "Band number" ), 1, QStringLiteral( "INPUT" ) ) );
57 addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "ZONES" ),
58 QObject::tr( "Zones layer" ) ) );
59 addParameter( new QgsProcessingParameterBand( QStringLiteral( "ZONES_BAND" ),
60 QObject::tr( "Zones band number" ), 1, QStringLiteral( "ZONES" ) ) );
61
62 std::unique_ptr< QgsProcessingParameterEnum > refParam = std::make_unique< QgsProcessingParameterEnum >( QStringLiteral( "REF_LAYER" ), QObject::tr( "Reference layer" ),
63 QStringList() << QObject::tr( "Input layer" ) << QObject::tr( "Zones layer" ), false, 0 );
64 refParam->setFlags( refParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
65 addParameter( refParam.release() );
66
67 addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT_TABLE" ),
68 QObject::tr( "Statistics" ), Qgis::ProcessingSourceType::Vector ) );
69
70 addOutput( new QgsProcessingOutputString( QStringLiteral( "EXTENT" ), QObject::tr( "Extent" ) ) );
71 addOutput( new QgsProcessingOutputString( QStringLiteral( "CRS_AUTHID" ), QObject::tr( "CRS authority identifier" ) ) );
72 addOutput( new QgsProcessingOutputNumber( QStringLiteral( "WIDTH_IN_PIXELS" ), QObject::tr( "Width in pixels" ) ) );
73 addOutput( new QgsProcessingOutputNumber( QStringLiteral( "HEIGHT_IN_PIXELS" ), QObject::tr( "Height in pixels" ) ) );
74 addOutput( new QgsProcessingOutputNumber( QStringLiteral( "TOTAL_PIXEL_COUNT" ), QObject::tr( "Total pixel count" ) ) );
75 addOutput( new QgsProcessingOutputNumber( QStringLiteral( "NODATA_PIXEL_COUNT" ), QObject::tr( "NoData pixel count" ) ) );
76}
77
78QString QgsRasterLayerZonalStatsAlgorithm::shortDescription() const
79{
80 return QObject::tr( "Calculates statistics for a raster layer's values, categorized by zones defined in another raster layer." );
81}
82
83QString QgsRasterLayerZonalStatsAlgorithm::shortHelpString() const
84{
85 return QObject::tr( "This algorithm calculates statistics for a raster layer's values, categorized by zones defined in another raster layer.\n\n"
86 "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"
87 "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"
88 "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." );
89}
90
91QgsRasterLayerZonalStatsAlgorithm *QgsRasterLayerZonalStatsAlgorithm::createInstance() const
92{
93 return new QgsRasterLayerZonalStatsAlgorithm();
94}
95
96bool QgsRasterLayerZonalStatsAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
97{
98 mRefLayer = static_cast< RefLayer >( parameterAsEnum( parameters, QStringLiteral( "REF_LAYER" ), context ) );
99
100 QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT" ), context );
101 const int band = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
102
103 if ( !layer )
104 throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT" ) ) );
105
106 mBand = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
107 if ( mBand < 1 || mBand > layer->bandCount() )
108 throw QgsProcessingException( QObject::tr( "Invalid band number for BAND (%1): Valid values for input raster are 1 to %2" ).arg( mBand )
109 .arg( layer->bandCount() ) );
110
111 mHasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( band );
112
113 QgsRasterLayer *zonesLayer = parameterAsRasterLayer( parameters, QStringLiteral( "ZONES" ), context );
114
115 if ( !zonesLayer )
116 throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "ZONES" ) ) );
117
118 mZonesBand = parameterAsInt( parameters, QStringLiteral( "ZONES_BAND" ), context );
119 if ( mZonesBand < 1 || mZonesBand > zonesLayer->bandCount() )
120 throw QgsProcessingException( QObject::tr( "Invalid band number for ZONES_BAND (%1): Valid values for input raster are 1 to %2" ).arg( mZonesBand )
121 .arg( zonesLayer->bandCount() ) );
122 mZonesHasNoDataValue = zonesLayer->dataProvider()->sourceHasNoDataValue( band );
123
124 mSourceDataProvider.reset( layer->dataProvider()->clone() );
125 mSourceInterface = mSourceDataProvider.get();
126 mZonesDataProvider.reset( zonesLayer->dataProvider()->clone() );
127 mZonesInterface = mZonesDataProvider.get();
128
129 switch ( mRefLayer )
130 {
131 case Source:
132 mCrs = layer->crs();
133 mRasterUnitsPerPixelX = layer->rasterUnitsPerPixelX();
134 mRasterUnitsPerPixelY = layer->rasterUnitsPerPixelY();
135 mLayerWidth = layer->width();
136 mLayerHeight = layer->height();
137 mExtent = layer->extent();
138
139 // add projector if necessary
140 if ( layer->crs() != zonesLayer->crs() )
141 {
142 mProjector = std::make_unique< QgsRasterProjector >();
143 mProjector->setInput( mZonesDataProvider.get() );
144 mProjector->setCrs( zonesLayer->crs(), layer->crs(), context.transformContext() );
145 mZonesInterface = mProjector.get();
146 }
147 break;
148
149 case Zones:
150 mCrs = zonesLayer->crs();
151 mRasterUnitsPerPixelX = zonesLayer->rasterUnitsPerPixelX();
152 mRasterUnitsPerPixelY = zonesLayer->rasterUnitsPerPixelY();
153 mLayerWidth = zonesLayer->width();
154 mLayerHeight = zonesLayer->height();
155 mExtent = zonesLayer->extent();
156
157 // add projector if necessary
158 if ( layer->crs() != zonesLayer->crs() )
159 {
160 mProjector = std::make_unique< QgsRasterProjector >();
161 mProjector->setInput( mSourceDataProvider.get() );
162 mProjector->setCrs( layer->crs(), zonesLayer->crs(), context.transformContext() );
163 mSourceInterface = mProjector.get();
164 }
165 break;
166 }
167
168 return true;
169}
170
171QVariantMap QgsRasterLayerZonalStatsAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
172{
173 QString areaUnit = QgsUnitTypes::toAbbreviatedString( QgsUnitTypes::distanceToAreaUnit( mCrs.mapUnits() ) );
174
175 QString tableDest;
176 std::unique_ptr< QgsFeatureSink > sink;
177 if ( parameters.contains( QStringLiteral( "OUTPUT_TABLE" ) ) && parameters.value( QStringLiteral( "OUTPUT_TABLE" ) ).isValid() )
178 {
179 QgsFields outFields;
180 outFields.append( QgsField( QStringLiteral( "zone" ), QMetaType::Type::Double, QString(), 20, 8 ) );
181 outFields.append( QgsField( areaUnit.isEmpty() ? "area" : areaUnit.replace( QStringLiteral( "²" ), QStringLiteral( "2" ) ), QMetaType::Type::Double, QString(), 20, 8 ) );
182 outFields.append( QgsField( QStringLiteral( "sum" ), QMetaType::Type::Double, QString(), 20, 8 ) );
183 outFields.append( QgsField( QStringLiteral( "count" ), QMetaType::Type::LongLong, QString(), 20 ) );
184 outFields.append( QgsField( QStringLiteral( "min" ), QMetaType::Type::Double, QString(), 20, 8 ) );
185 outFields.append( QgsField( QStringLiteral( "max" ), QMetaType::Type::Double, QString(), 20, 8 ) );
186 outFields.append( QgsField( QStringLiteral( "mean" ), QMetaType::Type::Double, QString(), 20, 8 ) );
187
188 sink.reset( parameterAsSink( parameters, QStringLiteral( "OUTPUT_TABLE" ), context, tableDest, outFields, Qgis::WkbType::NoGeometry, QgsCoordinateReferenceSystem() ) );
189 if ( !sink )
190 throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT_TABLE" ) ) );
191 }
192
193 struct StatCalculator
194 {
195 // only calculate cheap stats-- we cannot calculate stats which require holding values in memory -- because otherwise we'll end
196 // up trying to store EVERY pixel value from the input in memory
198 };
199 QHash<double, StatCalculator > zoneStats;
200 qgssize noDataCount = 0;
201
202 const qgssize layerSize = static_cast< qgssize >( mLayerWidth ) * static_cast< qgssize >( mLayerHeight );
205 const int nbBlocksWidth = static_cast< int>( std::ceil( 1.0 * mLayerWidth / maxWidth ) );
206 const int nbBlocksHeight = static_cast< int >( std::ceil( 1.0 * mLayerHeight / maxHeight ) );
207 const int nbBlocks = nbBlocksWidth * nbBlocksHeight;
208
209 QgsRasterIterator iter = mRefLayer == Source ? QgsRasterIterator( mSourceInterface )
210 : QgsRasterIterator( mZonesInterface );
211 iter.startRasterRead( mRefLayer == Source ? mBand : mZonesBand, mLayerWidth, mLayerHeight, mExtent );
212
213 int iterLeft = 0;
214 int iterTop = 0;
215 int iterCols = 0;
216 int iterRows = 0;
217 QgsRectangle blockExtent;
218 std::unique_ptr< QgsRasterBlock > rasterBlock;
219 std::unique_ptr< QgsRasterBlock > zonesRasterBlock;
220 bool isNoData = false;
221 while ( true )
222 {
223 if ( mRefLayer == Source )
224 {
225 if ( !iter.readNextRasterPart( mBand, iterCols, iterRows, rasterBlock, iterLeft, iterTop, &blockExtent ) )
226 break;
227
228 zonesRasterBlock.reset( mZonesInterface->block( mZonesBand, blockExtent, iterCols, iterRows ) );
229 }
230 else
231 {
232 if ( !iter.readNextRasterPart( mZonesBand, iterCols, iterRows, zonesRasterBlock, iterLeft, iterTop, &blockExtent ) )
233 break;
234
235 rasterBlock.reset( mSourceInterface->block( mBand, blockExtent, iterCols, iterRows ) );
236 }
237 if ( !zonesRasterBlock || !rasterBlock )
238 continue;
239
240 feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
241 if ( !rasterBlock->isValid() || rasterBlock->isEmpty() || !zonesRasterBlock->isValid() || zonesRasterBlock->isEmpty() )
242 continue;
243
244 for ( int row = 0; row < iterRows; row++ )
245 {
246 if ( feedback->isCanceled() )
247 break;
248
249 for ( int column = 0; column < iterCols; column++ )
250 {
251 const double value = rasterBlock->valueAndNoData( row, column, isNoData );
252 if ( mHasNoDataValue && isNoData )
253 {
254 noDataCount += 1;
255 continue;
256 }
257 const double zone = zonesRasterBlock->valueAndNoData( row, column, isNoData );
258 if ( mZonesHasNoDataValue && isNoData )
259 {
260 noDataCount += 1;
261 continue;
262 }
263 zoneStats[ zone ].s.addValue( value );
264 }
265 }
266 }
267
268 QVariantMap outputs;
269 outputs.insert( QStringLiteral( "EXTENT" ), mExtent.toString() );
270 outputs.insert( QStringLiteral( "CRS_AUTHID" ), mCrs.authid() );
271 outputs.insert( QStringLiteral( "WIDTH_IN_PIXELS" ), mLayerWidth );
272 outputs.insert( QStringLiteral( "HEIGHT_IN_PIXELS" ), mLayerHeight );
273 outputs.insert( QStringLiteral( "TOTAL_PIXEL_COUNT" ), layerSize );
274 outputs.insert( QStringLiteral( "NODATA_PIXEL_COUNT" ), noDataCount );
275
276 const double pixelArea = mRasterUnitsPerPixelX * mRasterUnitsPerPixelY;
277
278 for ( auto it = zoneStats.begin(); it != zoneStats.end(); ++it )
279 {
280 QgsFeature f;
281 it->s.finalize();
282 f.setAttributes( QgsAttributes() << it.key() << it->s.count() * pixelArea << it->s.sum() << it->s.count() <<
283 it->s.min() << it->s.max() << it->s.mean() );
284 if ( !sink->addFeature( f, QgsFeatureSink::FastInsert ) )
285 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT_TABLE" ) ) );
286 }
287 outputs.insert( QStringLiteral( "OUTPUT_TABLE" ), tableDest );
288
289 return outputs;
290}
291
292
294
295
296
@ 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:69
virtual QgsRectangle extent() const
Returns the extent of the layer.
QgsCoordinateReferenceSystem crs
Definition qgsmaplayer.h:82
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:6013