QGIS API Documentation  3.2.0-Bonn (bc43194)
qgsalgorithmrasterlayeruniquevalues.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsalgorithmrasterlayeruniquevalues.cpp
3  ---------------------
4  begin : April 2017
5  copyright : (C) 2017 by Mathieu Pellerin
6  email : nirvn dot asia 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 
22 
23 QString QgsRasterLayerUniqueValuesReportAlgorithm::name() const
24 {
25  return QStringLiteral( "rasterlayeruniquevaluesreport" );
26 }
27 
28 QString QgsRasterLayerUniqueValuesReportAlgorithm::displayName() const
29 {
30  return QObject::tr( "Raster layer unique values report" );
31 }
32 
33 QStringList QgsRasterLayerUniqueValuesReportAlgorithm::tags() const
34 {
35  return QObject::tr( "count,area,statistics" ).split( ',' );
36 }
37 
38 QString QgsRasterLayerUniqueValuesReportAlgorithm::group() const
39 {
40  return QObject::tr( "Raster analysis" );
41 }
42 
43 QString QgsRasterLayerUniqueValuesReportAlgorithm::groupId() const
44 {
45  return QStringLiteral( "rasteranalysis" );
46 }
47 
48 void QgsRasterLayerUniqueValuesReportAlgorithm::initAlgorithm( const QVariantMap & )
49 {
50  addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT" ),
51  QObject::tr( "Input layer" ) ) );
52  addParameter( new QgsProcessingParameterBand( QStringLiteral( "BAND" ),
53  QObject::tr( "Band number" ), 1, QStringLiteral( "INPUT" ) ) );
54  addParameter( new QgsProcessingParameterFileDestination( QStringLiteral( "OUTPUT_HTML_FILE" ),
55  QObject::tr( "Unique values report" ), QObject::tr( "HTML files (*.html)" ), QVariant(), true ) );
56  addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT_TABLE" ),
57  QObject::tr( "Unique values table" ), QgsProcessing::TypeVector, QVariant(), true, false ) );
58 
59  addOutput( new QgsProcessingOutputString( QStringLiteral( "EXTENT" ), QObject::tr( "Extent" ) ) );
60  addOutput( new QgsProcessingOutputString( QStringLiteral( "CRS_AUTHID" ), QObject::tr( "CRS authority identifier" ) ) );
61  addOutput( new QgsProcessingOutputNumber( QStringLiteral( "WIDTH_IN_PIXELS" ), QObject::tr( "Width in pixels" ) ) );
62  addOutput( new QgsProcessingOutputNumber( QStringLiteral( "HEIGHT_IN_PIXELS" ), QObject::tr( "Height in pixels" ) ) );
63  addOutput( new QgsProcessingOutputNumber( QStringLiteral( "TOTAL_PIXEL_COUNT" ), QObject::tr( "Total pixel count" ) ) );
64  addOutput( new QgsProcessingOutputNumber( QStringLiteral( "NODATA_PIXEL_COUNT" ), QObject::tr( "NODATA pixel count" ) ) );
65 }
66 
67 QString QgsRasterLayerUniqueValuesReportAlgorithm::shortHelpString() const
68 {
69  return QObject::tr( "This algorithm returns the count and area of each unique value in a given raster layer." );
70 }
71 
72 QgsRasterLayerUniqueValuesReportAlgorithm *QgsRasterLayerUniqueValuesReportAlgorithm::createInstance() const
73 {
74  return new QgsRasterLayerUniqueValuesReportAlgorithm();
75 }
76 
77 bool QgsRasterLayerUniqueValuesReportAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
78 {
79  QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT" ), context );
80  int band = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
81 
82  if ( !layer )
83  throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT" ) ) );
84 
85  mBand = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
86  if ( mBand < 1 || mBand > layer->bandCount() )
87  throw QgsProcessingException( QObject::tr( "Invalid band number for BAND (%1): Valid values for input raster are 1 to %2" ).arg( mBand )
88  .arg( layer->bandCount() ) );
89 
90  mInterface.reset( layer->dataProvider()->clone() );
91  mHasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( band );
92  mLayerWidth = layer->width();
93  mLayerHeight = layer->height();
94  mExtent = layer->extent();
95  mCrs = layer->crs();
96  mRasterUnitsPerPixelX = layer->rasterUnitsPerPixelX();
97  mRasterUnitsPerPixelY = layer->rasterUnitsPerPixelY();
98  mSource = layer->source();
99 
100  return true;
101 }
102 
103 QVariantMap QgsRasterLayerUniqueValuesReportAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
104 {
105  QString outputFile = parameterAsFileOutput( parameters, QStringLiteral( "OUTPUT_HTML_FILE" ), context );
106 
107  QString areaUnit = QgsUnitTypes::toAbbreviatedString( QgsUnitTypes::distanceToAreaUnit( mCrs.mapUnits() ) );
108 
109  QString tableDest;
110  std::unique_ptr< QgsFeatureSink > sink;
111  if ( parameters.contains( QStringLiteral( "OUTPUT_TABLE" ) ) && parameters.value( QStringLiteral( "OUTPUT_TABLE" ) ).isValid() )
112  {
113  QgsFields outFields;
114  outFields.append( QgsField( QStringLiteral( "value" ), QVariant::Double, QString(), 20, 8 ) );
115  outFields.append( QgsField( QStringLiteral( "count" ), QVariant::Int, QString(), 20 ) );
116  outFields.append( QgsField( areaUnit.replace( QStringLiteral( "²" ), QStringLiteral( "2" ) ), QVariant::Double, QString(), 20, 8 ) );
117  sink.reset( parameterAsSink( parameters, QStringLiteral( "OUTPUT_TABLE" ), context, tableDest, outFields, QgsWkbTypes::NoGeometry, QgsCoordinateReferenceSystem() ) );
118  if ( !sink )
119  throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT_TABLE" ) ) );
120  }
121 
122  QHash< double, qgssize > uniqueValues;
123  qgssize noDataCount = 0;
124 
125  qgssize layerSize = static_cast< qgssize >( mLayerWidth ) * static_cast< qgssize >( mLayerHeight );
128  int nbBlocksWidth = std::ceil( 1.0 * mLayerWidth / maxWidth );
129  int nbBlocksHeight = std::ceil( 1.0 * mLayerHeight / maxHeight );
130  int nbBlocks = nbBlocksWidth * nbBlocksHeight;
131 
132  QgsRasterIterator iter( mInterface.get() );
133  iter.startRasterRead( mBand, mLayerWidth, mLayerHeight, mExtent );
134 
135  int iterLeft = 0;
136  int iterTop = 0;
137  int iterCols = 0;
138  int iterRows = 0;
139  std::unique_ptr< QgsRasterBlock > rasterBlock;
140  while ( iter.readNextRasterPart( mBand, iterCols, iterRows, rasterBlock, iterLeft, iterTop ) )
141  {
142  feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
143  for ( int row = 0; row < iterRows; row++ )
144  {
145  if ( feedback->isCanceled() )
146  break;
147  for ( int column = 0; column < iterCols; column++ )
148  {
149  if ( mHasNoDataValue && rasterBlock->isNoData( row, column ) )
150  {
151  noDataCount += 1;
152  }
153  else
154  {
155  double value = rasterBlock->value( row, column );
156  uniqueValues[ value ]++;
157  }
158  }
159  }
160  }
161 
162  QMap< double, qgssize > sortedUniqueValues;
163  for ( auto it = uniqueValues.constBegin(); it != uniqueValues.constEnd(); ++it )
164  {
165  sortedUniqueValues.insert( it.key(), it.value() );
166  }
167 
168  QVariantMap outputs;
169  outputs.insert( QStringLiteral( "EXTENT" ), mExtent.toString() );
170  outputs.insert( QStringLiteral( "CRS_AUTHID" ), mCrs.authid() );
171  outputs.insert( QStringLiteral( "WIDTH_IN_PIXELS" ), mLayerWidth );
172  outputs.insert( QStringLiteral( "HEIGHT_IN_PIXELS" ), mLayerHeight );
173  outputs.insert( QStringLiteral( "TOTAL_PIXEL_COUNT" ), layerSize );
174  outputs.insert( QStringLiteral( "NODATA_PIXEL_COUNT" ), noDataCount );
175 
176  double pixelArea = mRasterUnitsPerPixelX * mRasterUnitsPerPixelY;
177 
178  if ( !outputFile.isEmpty() )
179  {
180  QFile file( outputFile );
181  if ( file.open( QIODevice::WriteOnly | QIODevice::Text ) )
182  {
183  const QString encodedAreaUnit = QgsStringUtils::ampersandEncode( areaUnit );
184 
185  QTextStream out( &file );
186  out << QStringLiteral( "<html><head><meta http-equiv=\"Content-Type\" content=\"text/html;charset=utf-8\"/></head><body>\n" );
187  out << QStringLiteral( "<p>%1: %2 (%3 %4)</p>\n" ).arg( QObject::tr( "Analyzed file" ), mSource, QObject::tr( "band" ) ).arg( mBand );
188  out << QObject::tr( "<p>%1: %2</p>\n" ).arg( QObject::tr( "Extent" ), mExtent.toString() );
189  out << QObject::tr( "<p>%1: %2 (%3)</p>\n" ).arg( QObject::tr( "Projection" ), mCrs.description(), mCrs.authid() );
190  out << QObject::tr( "<p>%1: %2 (%3 %4)</p>\n" ).arg( QObject::tr( "Width in pixels" ) ).arg( mLayerWidth ).arg( QObject::tr( "units per pixel" ) ).arg( mRasterUnitsPerPixelX );
191  out << QObject::tr( "<p>%1: %2 (%3 %4)</p>\n" ).arg( QObject::tr( "Height in pixels" ) ).arg( mLayerHeight ).arg( QObject::tr( "units per pixel" ) ).arg( mRasterUnitsPerPixelY );
192  out << QObject::tr( "<p>%1: %2</p>\n" ).arg( QObject::tr( "Total pixel count" ) ).arg( layerSize );
193  if ( mHasNoDataValue )
194  out << QObject::tr( "<p>%1: %2</p>\n" ).arg( QObject::tr( "NODATA pixel count" ) ).arg( noDataCount );
195  out << QStringLiteral( "<table><tr><td>%1</td><td>%2</td><td>%3 (%4)</td></tr>\n" ).arg( QObject::tr( "Value" ), QObject::tr( "Pixel count" ), QObject::tr( "Area" ), encodedAreaUnit );
196 
197  for ( auto it = sortedUniqueValues.constBegin(); it != sortedUniqueValues.constEnd(); ++it )
198  {
199  double area = it.value() * pixelArea;
200  out << QStringLiteral( "<tr><td>%1</td><td>%2</td><td>%3</td></tr>\n" ).arg( it.key() ).arg( it.value() ).arg( QString::number( area, 'g', 16 ) );
201  }
202  out << QStringLiteral( "</table>\n</body></html>" );
203  outputs.insert( QStringLiteral( "OUTPUT_HTML_FILE" ), outputFile );
204  }
205  }
206 
207  if ( sink )
208  {
209  for ( auto it = sortedUniqueValues.constBegin(); it != sortedUniqueValues.constEnd(); ++it )
210  {
211  QgsFeature f;
212  double area = it.value() * pixelArea;
213  f.setAttributes( QgsAttributes() << it.key() << it.value() << area );
214  sink->addFeature( f, QgsFeatureSink::FastInsert );
215  }
216  outputs.insert( QStringLiteral( "OUTPUT_TABLE" ), tableDest );
217  }
218 
219  return outputs;
220 }
221 
222 
224 
225 
226 
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...
Base class for providing feedback from a processing algorithm.
Iterator for sequentially processing raster cells.
int bandCount() const
Returns the number of bands in this layer.
This class provides qgis with the ability to render raster datasets onto the mapcanvas.
double rasterUnitsPerPixelX() const
Returns the number of raster units per each raster pixel in X axis.
QgsRasterInterface * clone() const override=0
Clone itself, create deep copy.
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
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:42
void setAttributes(const QgsAttributes &attrs)
Sets the feature&#39;s attributes.
Definition: qgsfeature.cpp:127
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:62
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 layer&#39;s 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.
static const int DEFAULT_MAXIMUM_TILE_HEIGHT
Default maximum tile height.
Custom exception class for processing related exceptions.
Definition: qgsexception.h:82
QgsCoordinateReferenceSystem crs() const
Returns the layer&#39;s spatial reference system.
bool append(const QgsField &field, FieldOrigin origin=OriginProvider, int originIndex=-1)
Append 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:48
A generic file based destination parameter, for specifying the destination path for a file (non-map l...
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:510
static QString ampersandEncode(const QString &string)
Makes a raw string safe for inclusion as a HTML/XML string literal.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
QString source() const
Returns the source for the layer.
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:58
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.
void startRasterRead(int bandNumber, int nCols, int nRows, const QgsRectangle &extent, QgsRasterBlockFeedback *feedback=nullptr)
Start reading of raster band.