QGIS API Documentation 3.37.0-Master (fdefdf9c27f)
qgsalgorithmcalculateoverlaps.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsalgorithmcalculateoverlaps.cpp
3 ------------------
4 begin : May 2019
5 copyright : (C) 2019 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 "qgsvectorlayer.h"
20#include "qgsgeometryengine.h"
21#include "qgsdistancearea.h"
22#include "qgsspatialindex.h"
23
25
26QString QgsCalculateVectorOverlapsAlgorithm::name() const
27{
28 return QStringLiteral( "calculatevectoroverlaps" );
29}
30
31QString QgsCalculateVectorOverlapsAlgorithm::displayName() const
32{
33 return QObject::tr( "Overlap analysis" );
34}
35
36QStringList QgsCalculateVectorOverlapsAlgorithm::tags() const
37{
38 return QObject::tr( "vector,overlay,area,percentage,intersection" ).split( ',' );
39}
40
41QString QgsCalculateVectorOverlapsAlgorithm::group() const
42{
43 return QObject::tr( "Vector analysis" );
44}
45
46QString QgsCalculateVectorOverlapsAlgorithm::groupId() const
47{
48 return QStringLiteral( "vectoranalysis" );
49}
50
51void QgsCalculateVectorOverlapsAlgorithm::initAlgorithm( const QVariantMap & )
52{
53 addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT" ), QObject::tr( "Input layer" ), QList< int >() << static_cast< int >( Qgis::ProcessingSourceType::VectorPolygon ) ) );
54 addParameter( new QgsProcessingParameterMultipleLayers( QStringLiteral( "LAYERS" ), QObject::tr( "Overlay layers" ), Qgis::ProcessingSourceType::VectorPolygon ) );
55 addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Overlap" ) ) );
56
57 std::unique_ptr< QgsProcessingParameterNumber > gridSize = std::make_unique< QgsProcessingParameterNumber >( QStringLiteral( "GRID_SIZE" ),
58 QObject::tr( "Grid size" ), Qgis::ProcessingNumberParameterType::Double, QVariant(), true, 0 );
59 gridSize->setFlags( gridSize->flags() | Qgis::ProcessingParameterFlag::Advanced );
60 addParameter( gridSize.release() );
61}
62
63QIcon QgsCalculateVectorOverlapsAlgorithm::icon() const
64{
65 return QgsApplication::getThemeIcon( QStringLiteral( "/algorithms/mAlgorithmClip.svg" ) );
66}
67
68QString QgsCalculateVectorOverlapsAlgorithm::svgIconPath() const
69{
70 return QgsApplication::iconPath( QStringLiteral( "/algorithms/mAlgorithmClip.svg" ) );
71}
72
73QString QgsCalculateVectorOverlapsAlgorithm::shortHelpString() const
74{
75 return QObject::tr( "This algorithm calculates the area and percentage cover by which features from an input layer "
76 "are overlapped by features from a selection of overlay layers.\n\n"
77 "New attributes are added to the output layer reporting the total area of overlap and percentage of the input feature overlapped "
78 "by each of the selected overlay layers." );
79}
80
81QgsCalculateVectorOverlapsAlgorithm *QgsCalculateVectorOverlapsAlgorithm::createInstance() const
82{
83 return new QgsCalculateVectorOverlapsAlgorithm();
84}
85
86bool QgsCalculateVectorOverlapsAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
87{
88 mSource.reset( parameterAsSource( parameters, QStringLiteral( "INPUT" ), context ) );
89 if ( !mSource )
90 throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT" ) ) );
91
92 mOutputFields = mSource->fields();
93
94 const QList< QgsMapLayer * > layers = parameterAsLayerList( parameters, QStringLiteral( "LAYERS" ), context );
95 mOverlayerSources.reserve( layers.size() );
96 mLayerNames.reserve( layers.size() );
97 for ( QgsMapLayer *layer : layers )
98 {
99 if ( QgsVectorLayer *vl = qobject_cast< QgsVectorLayer * >( layer ) )
100 {
101 mLayerNames << layer->name();
102 mOverlayerSources.emplace_back( std::make_unique< QgsVectorLayerFeatureSource >( vl ) );
103 mOutputFields.append( QgsField( QStringLiteral( "%1_area" ).arg( vl->name() ), QVariant::Double ) );
104 mOutputFields.append( QgsField( QStringLiteral( "%1_pc" ).arg( vl->name() ), QVariant::Double ) );
105 }
106 }
107
108 mOutputType = mSource->wkbType();
109 mCrs = mSource->sourceCrs();
110 mInputCount = mSource->featureCount();
111 mInputFeatures = mSource->getFeatures();
112 return true;
113}
114
115QVariantMap QgsCalculateVectorOverlapsAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
116{
117 QString destId;
118 std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, destId, mOutputFields,
119 mOutputType, mCrs ) );
120 if ( !sink )
121 throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );
122
123 // build a spatial index for each constraint layer for speed. We also store input constraint geometries here,
124 // to avoid refetching and projecting them later
125 QList< QgsSpatialIndex > spatialIndices;
126 spatialIndices.reserve( mLayerNames.size() );
127 auto nameIt = mLayerNames.constBegin();
128 for ( auto sourceIt = mOverlayerSources.begin(); sourceIt != mOverlayerSources.end(); ++sourceIt, ++nameIt )
129 {
130 feedback->pushInfo( QObject::tr( "Preparing %1" ).arg( *nameIt ) );
131 const QgsFeatureIterator featureIt = ( *sourceIt )->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( QgsAttributeList() ).setDestinationCrs( mCrs, context.transformContext() ).setInvalidGeometryCheck( context.invalidGeometryCheck() ).setInvalidGeometryCallback( context.invalidGeometryCallback() ) );
132 spatialIndices << QgsSpatialIndex( featureIt, feedback, QgsSpatialIndex::FlagStoreFeatureGeometries );
133 }
134
136 da.setSourceCrs( mCrs, context.transformContext() );
137 da.setEllipsoid( context.ellipsoid() );
138
139 QgsGeometryParameters geometryParameters;
140 if ( parameters.value( QStringLiteral( "GRID_SIZE" ) ).isValid() )
141 {
142 geometryParameters.setGridSize( parameterAsDouble( parameters, QStringLiteral( "GRID_SIZE" ), context ) );
143 }
144
145 // loop through input
146 const double step = mInputCount > 0 ? 100.0 / mInputCount : 0;
147 long i = 0;
148 QgsFeature feature;
149 while ( mInputFeatures.nextFeature( feature ) )
150 {
151 if ( feedback->isCanceled() )
152 break;
153
154 QgsAttributes outAttributes = feature.attributes();
155 if ( feature.hasGeometry() && !qgsDoubleNear( feature.geometry().area(), 0.0 ) )
156 {
157 const QgsGeometry inputGeom = feature.geometry();
158 const double inputArea = da.measureArea( inputGeom );
159
160 // prepare for lots of intersection tests (for speed)
161 std::unique_ptr< QgsGeometryEngine > bufferGeomEngine( QgsGeometry::createGeometryEngine( inputGeom.constGet() ) );
162 bufferGeomEngine->prepareGeometry();
163
164 // calculate overlap attributes
165 auto spatialIteratorIt = spatialIndices.begin();
166 for ( auto it = mOverlayerSources.begin(); it != mOverlayerSources.end(); ++ it, ++spatialIteratorIt )
167 {
168 if ( feedback->isCanceled() )
169 break;
170
171 const QgsSpatialIndex &index = *spatialIteratorIt;
172 const QList<QgsFeatureId> matches = index.intersects( inputGeom.boundingBox() );
173 QVector< QgsGeometry > intersectingGeoms;
174 intersectingGeoms.reserve( matches.count() );
175 for ( const QgsFeatureId match : matches )
176 {
177 if ( feedback->isCanceled() )
178 break;
179
180 const QgsGeometry overlayGeometry = index.geometry( match );
181 if ( bufferGeomEngine->intersects( overlayGeometry.constGet() ) )
182 {
183 intersectingGeoms.append( overlayGeometry );
184 }
185 }
186
187 if ( feedback->isCanceled() )
188 break;
189
190 // dissolve intersecting features, calculate total area of them within our buffer
191 const QgsGeometry overlayDissolved = QgsGeometry::unaryUnion( intersectingGeoms, geometryParameters );
192
193 if ( feedback->isCanceled() )
194 break;
195
196 const QgsGeometry overlayIntersection = inputGeom.intersection( overlayDissolved, geometryParameters );
197
198 const double overlayArea = da.measureArea( overlayIntersection );
199 outAttributes.append( overlayArea );
200 outAttributes.append( 100 * overlayArea / inputArea );
201 }
202 }
203 else
204 {
205 // input feature has no geometry
206 for ( auto it = mOverlayerSources.begin(); it != mOverlayerSources.end(); ++ it )
207 {
208 outAttributes.append( QVariant() );
209 outAttributes.append( QVariant() );
210 }
211 }
212
213 feature.setAttributes( outAttributes );
214 if ( !sink->addFeature( feature, QgsFeatureSink::FastInsert ) )
215 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
216
217 i++;
218 feedback->setProgress( i * step );
219 }
220
221 QVariantMap outputs;
222 outputs.insert( QStringLiteral( "OUTPUT" ), destId );
223 return outputs;
224}
225
@ VectorPolygon
Vector polygon layers.
@ Advanced
Parameter is an advanced parameter which should be hidden from users by default.
static QIcon getThemeIcon(const QString &name, const QColor &fillColor=QColor(), const QColor &strokeColor=QColor())
Helper to get a theme icon.
static QString iconPath(const QString &iconFile)
Returns path to the desired icon file.
A vector of attributes.
Definition: qgsattributes.h:59
A general purpose distance and area calculator, capable of performing ellipsoid based calculations.
double measureArea(const QgsGeometry &geometry) const
Measures the area of a geometry.
void setSourceCrs(const QgsCoordinateReferenceSystem &crs, const QgsCoordinateTransformContext &context)
Sets source spatial reference system crs.
bool setEllipsoid(const QString &ellipsoid)
Sets the ellipsoid by its acronym.
Wrapper for iterator of features from vector data provider or vector layer.
This class wraps a request for features to a vector layer (or directly its vector data provider).
QgsFeatureRequest & setInvalidGeometryCheck(Qgis::InvalidGeometryCheck check)
Sets invalid geometry checking behavior.
QgsFeatureRequest & setInvalidGeometryCallback(const std::function< void(const QgsFeature &)> &callback)
Sets a callback function to use when encountering an invalid geometry and invalidGeometryCheck() is s...
@ 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:56
QgsAttributes attributes
Definition: qgsfeature.h:65
void setAttributes(const QgsAttributes &attrs)
Sets the feature's attributes.
Definition: qgsfeature.cpp:160
QgsGeometry geometry
Definition: qgsfeature.h:67
bool hasGeometry() const
Returns true if the feature has an associated geometry.
Definition: qgsfeature.cpp:230
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
Encapsulates parameters under which a geometry operation is performed.
Definition: qgsgeometry.h:109
void setGridSize(double size)
Sets the grid size which will be used to snap vertices of a geometry.
Definition: qgsgeometry.h:134
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:162
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
double area() const
Returns the planar, 2-dimensional area of the geometry.
QgsGeometry intersection(const QgsGeometry &geometry, const QgsGeometryParameters &parameters=QgsGeometryParameters()) const
Returns a geometry representing the points shared by this geometry and other.
static QgsGeometry unaryUnion(const QVector< QgsGeometry > &geometries, const QgsGeometryParameters &parameters=QgsGeometryParameters())
Compute the unary union on a list of geometries.
QgsRectangle boundingBox() const
Returns the bounding box of the geometry.
static QgsGeometryEngine * createGeometryEngine(const QgsAbstractGeometry *geometry, double precision=0.0)
Creates and returns a new geometry engine representing the specified geometry using precision on a gr...
Base class for all map layer types.
Definition: qgsmaplayer.h:75
Contains information about the context in which a processing algorithm is executed.
QgsCoordinateTransformContext transformContext() const
Returns the coordinate transform context.
std::function< void(const QgsFeature &) > invalidGeometryCallback(QgsFeatureSource *source=nullptr) const
Returns the callback function to use when encountering an invalid geometry and invalidGeometryCheck()...
QString ellipsoid() const
Returns the ellipsoid to use for distance and area calculations.
Qgis::InvalidGeometryCheck invalidGeometryCheck() const
Returns the behavior used for checking invalid geometries in input layers.
Custom exception class for processing related exceptions.
Definition: qgsexception.h:83
Base class for providing feedback from a processing algorithm.
virtual void pushInfo(const QString &info)
Pushes a general informational message from the algorithm.
A feature sink output for processing algorithms.
An input feature source (such as vector layers) parameter for processing algorithms.
A parameter for processing algorithms which accepts multiple map layers.
A spatial index for QgsFeature objects.
@ FlagStoreFeatureGeometries
Indicates that the spatial index should also store feature geometries. This requires more memory,...
QList< QgsFeatureId > intersects(const QgsRectangle &rectangle) const
Returns a list of features with a bounding box which intersects the specified rectangle.
QgsGeometry geometry(QgsFeatureId id) const
Returns the stored geometry for the indexed feature with matching id.
Represents a vector layer which manages a vector based data sets.
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:5207
qint64 QgsFeatureId
64 bit feature ids negative numbers are used for uncommitted/newly added features
Definition: qgsfeatureid.h:28
QList< int > QgsAttributeList
Definition: qgsfield.h:27