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