QGIS API Documentation  3.10.0-A Coruña (6c816b4204)
qgsalgorithmzonalhistogram.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsalgorithmzonalhistogram.cpp
3  ---------------------
4  begin : May, 2018
5  copyright : (C) 2018 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 "qgsrasteranalysisutils.h"
20 #include "qgslogger.h"
21 
23 
24 QString QgsZonalHistogramAlgorithm::name() const
25 {
26  return QStringLiteral( "zonalhistogram" );
27 }
28 
29 QString QgsZonalHistogramAlgorithm::displayName() const
30 {
31  return QObject::tr( "Zonal histogram" );
32 }
33 
34 QStringList QgsZonalHistogramAlgorithm::tags() const
35 {
36  return QObject::tr( "raster,unique,values,count,area,statistics" ).split( ',' );
37 }
38 
39 QString QgsZonalHistogramAlgorithm::group() const
40 {
41  return QObject::tr( "Raster analysis" );
42 }
43 
44 QString QgsZonalHistogramAlgorithm::groupId() const
45 {
46  return QStringLiteral( "rasteranalysis" );
47 }
48 
49 void QgsZonalHistogramAlgorithm::initAlgorithm( const QVariantMap & )
50 {
51  addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT_RASTER" ),
52  QObject::tr( "Raster layer" ) ) );
53  addParameter( new QgsProcessingParameterBand( QStringLiteral( "RASTER_BAND" ),
54  QObject::tr( "Band number" ), 1, QStringLiteral( "INPUT_RASTER" ) ) );
55 
56  addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT_VECTOR" ),
57  QObject::tr( "Vector layer containing zones" ), QList< int >() << QgsProcessing::TypeVectorPolygon ) );
58 
59  addParameter( new QgsProcessingParameterString( QStringLiteral( "COLUMN_PREFIX" ), QObject::tr( "Output column prefix" ), QStringLiteral( "HISTO_" ), false, true ) );
60 
61  addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Output zones" ), QgsProcessing::TypeVectorPolygon ) );
62 }
63 
64 QString QgsZonalHistogramAlgorithm::shortHelpString() const
65 {
66  return QObject::tr( "This algorithm appends fields representing counts of each unique value from a raster layer contained within zones defined as polygons." );
67 }
68 
69 QgsZonalHistogramAlgorithm *QgsZonalHistogramAlgorithm::createInstance() const
70 {
71  return new QgsZonalHistogramAlgorithm();
72 }
73 
74 bool QgsZonalHistogramAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
75 {
76  QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT_RASTER" ), context );
77  if ( !layer )
78  throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "INPUT_RASTER" ) ) );
79 
80  mRasterBand = parameterAsInt( parameters, QStringLiteral( "RASTER_BAND" ), context );
81  mHasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( mRasterBand );
82  mNodataValue = layer->dataProvider()->sourceNoDataValue( mRasterBand );
83  mRasterInterface.reset( layer->dataProvider()->clone() );
84  mRasterExtent = layer->extent();
85  mCrs = layer->crs();
86  mCellSizeX = std::abs( layer->rasterUnitsPerPixelX() );
87  mCellSizeY = std::abs( layer->rasterUnitsPerPixelX() );
88  mNbCellsXProvider = mRasterInterface->xSize();
89  mNbCellsYProvider = mRasterInterface->ySize();
90 
91  return true;
92 }
93 
94 QVariantMap QgsZonalHistogramAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
95 {
96 
97  std::unique_ptr< QgsFeatureSource > zones( parameterAsSource( parameters, QStringLiteral( "INPUT_VECTOR" ), context ) );
98  if ( !zones )
99  throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT_VECTOR" ) ) );
100 
101  long count = zones->featureCount();
102  double step = count > 0 ? 100.0 / count : 1;
103  long current = 0;
104 
105  QList< double > uniqueValues;
106  QMap< QgsFeatureId, QHash< double, qgssize > > featuresUniqueValues;
107 
108  // First loop through the zones to build up a list of unique values across all zones to determine sink fields list
109  QgsFeatureRequest request;
110  request.setNoAttributes();
111  if ( zones->sourceCrs() != mCrs )
112  {
113  request.setDestinationCrs( mCrs, context.transformContext() );
114  }
115  QgsFeatureIterator it = zones->getFeatures( request );
116  QgsFeature f;
117  while ( it.nextFeature( f ) )
118  {
119  if ( feedback && feedback->isCanceled() )
120  {
121  break;
122  }
123  feedback->setProgress( current * step );
124 
125  if ( !f.hasGeometry() )
126  {
127  current++;
128  continue;
129  }
130 
131  QgsGeometry featureGeometry = f.geometry();
132  QgsRectangle featureRect = featureGeometry.boundingBox().intersect( mRasterExtent );
133  if ( featureRect.isEmpty() )
134  {
135  current++;
136  continue;
137  }
138 
139  int nCellsX, nCellsY;
140  QgsRectangle rasterBlockExtent;
141  QgsRasterAnalysisUtils::cellInfoForBBox( mRasterExtent, featureRect, mCellSizeX, mCellSizeY, nCellsX, nCellsY, mNbCellsXProvider, mNbCellsYProvider, rasterBlockExtent );
142 
143  QHash< double, qgssize > fUniqueValues;
144  QgsRasterAnalysisUtils::statisticsFromMiddlePointTest( mRasterInterface.get(), mRasterBand, featureGeometry, nCellsX, nCellsY, mCellSizeX, mCellSizeY,
145  rasterBlockExtent, [ &fUniqueValues]( double value ) { fUniqueValues[value]++; }, false );
146 
147  if ( fUniqueValues.count() < 1 )
148  {
149  // The cell resolution is probably larger than the polygon area. We switch to slower precise pixel - polygon intersection in this case
150  // TODO: eventually deal with weight if needed
151  QgsRasterAnalysisUtils::statisticsFromPreciseIntersection( mRasterInterface.get(), mRasterBand, featureGeometry, nCellsX, nCellsY, mCellSizeX, mCellSizeY,
152  rasterBlockExtent, [ &fUniqueValues]( double value, double ) { fUniqueValues[value]++; }, false );
153  }
154 
155  for ( auto it = fUniqueValues.constBegin(); it != fUniqueValues.constEnd(); ++it )
156  {
157  if ( uniqueValues.indexOf( it.key() ) == -1 )
158  {
159  uniqueValues << it.key();
160  }
161  featuresUniqueValues[f.id()][it.key()] += it.value();
162  }
163 
164  current++;
165  }
166 
167  std::sort( uniqueValues.begin(), uniqueValues.end() );
168 
169  QString fieldPrefix = parameterAsString( parameters, QStringLiteral( "COLUMN_PREFIX" ), context );
170  QgsFields newFields;
171  for ( auto it = uniqueValues.constBegin(); it != uniqueValues.constEnd(); ++it )
172  {
173  newFields.append( QgsField( QStringLiteral( "%1%2" ).arg( fieldPrefix, mHasNoDataValue && *it == mNodataValue ? QStringLiteral( "NODATA" ) : QString::number( *it ) ), QVariant::LongLong, QString(), -1, 0 ) );
174  }
175  QgsFields fields = QgsProcessingUtils::combineFields( zones->fields(), newFields );
176 
177  QString dest;
178  std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, fields,
179  zones->wkbType(), zones->sourceCrs() ) );
180  if ( !sink )
181  throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );
182 
183  it = zones->getFeatures( QgsFeatureRequest() );
184  while ( it.nextFeature( f ) )
185  {
186  QgsAttributes attributes = f.attributes();
187  QHash< double, qgssize > fUniqueValues = featuresUniqueValues.value( f.id() );
188  for ( auto it = uniqueValues.constBegin(); it != uniqueValues.constEnd(); ++it )
189  {
190  attributes += fUniqueValues.value( *it, 0 );
191  }
192 
193  QgsFeature outputFeature;
194  outputFeature.setGeometry( f.geometry() );
195  outputFeature.setAttributes( attributes );
196 
197  sink->addFeature( outputFeature, QgsFeatureSink::FastInsert );
198  }
199 
200  QVariantMap outputs;
201  outputs.insert( QStringLiteral( "OUTPUT" ), dest );
202  return outputs;
203 }
204 
206 
207 
208 
QgsFeatureRequest & setDestinationCrs(const QgsCoordinateReferenceSystem &crs, const QgsCoordinateTransformContext &context)
Sets the destination crs for feature&#39;s geometries.
QgsFeatureId id
Definition: qgsfeature.h:64
Wrapper for iterator of features from vector data provider or vector layer.
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.
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
Container of fields for a vector layer.
Definition: qgsfields.h:42
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:122
void setAttributes(const QgsAttributes &attrs)
Sets the feature&#39;s attributes.
Definition: qgsfeature.cpp:127
A raster band parameter for Processing algorithms.
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
Definition: qgsfeature.h:55
bool hasGeometry() const
Returns true if the feature has an associated geometry.
Definition: qgsfeature.cpp:197
virtual double sourceNoDataValue(int bandNo) const
Value representing no data value.
A feature sink output for processing algorithms.
virtual QgsRectangle extent() const
Returns the extent of the layer.
QgsRasterDataProvider * dataProvider() override
Returns the source data provider.
A raster layer parameter for processing algorithms.
bool isEmpty() const
Returns true if the rectangle is empty.
Definition: qgsrectangle.h:426
QgsFeatureRequest & setNoAttributes()
Set that no attributes will be fetched.
QgsRectangle intersect(const QgsRectangle &rect) const
Returns the intersection with the given rectangle.
Definition: qgsrectangle.h:312
This class wraps a request for features to a vector layer (or directly its vector data provider)...
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
Vector polygon layers.
Definition: qgsprocessing.h:50
Encapsulate a field in an attribute table or data source.
Definition: qgsfield.h:48
virtual bool sourceHasNoDataValue(int bandNo) const
Returns true if source band has no data value.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
An input feature source (such as vector layers) parameter for processing algorithms.
QgsRectangle boundingBox() const
Returns the bounding box of the geometry.
void setGeometry(const QgsGeometry &geometry)
Set the feature&#39;s geometry.
Definition: qgsfeature.cpp:137
QgsGeometry geometry
Definition: qgsfeature.h:67
bool nextFeature(QgsFeature &f)
static QgsFields combineFields(const QgsFields &fieldsA, const QgsFields &fieldsB, const QString &fieldsBPrefix=QString())
Combines two field lists, avoiding duplicate field names (in a case-insensitive manner).
A vector of attributes.
Definition: qgsattributes.h:57
Contains information about the context in which a processing algorithm is executed.
A string parameter for processing algorithms.
QgsAttributes attributes
Definition: qgsfeature.h:65
QgsCoordinateReferenceSystem crs
Definition: qgsmaplayer.h:86