QGIS API Documentation 3.40.0-Bratislava (b56115d8743)
Loading...
Searching...
No Matches
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() ), QMetaType::Type::Double ) );
104 mOutputFields.append( QgsField( QStringLiteral( "%1_pc" ).arg( vl->name() ), QMetaType::Type::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 ) );
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
159 double inputArea = 0;
160 try
161 {
162 inputArea = da.measureArea( inputGeom );
163 }
164 catch ( QgsCsException & )
165 {
166 throw QgsProcessingException( QObject::tr( "An error occurred while calculating feature area" ) );
167 }
168
169 // prepare for lots of intersection tests (for speed)
170 std::unique_ptr< QgsGeometryEngine > bufferGeomEngine( QgsGeometry::createGeometryEngine( inputGeom.constGet() ) );
171 bufferGeomEngine->prepareGeometry();
172
173 // calculate overlap attributes
174 auto spatialIteratorIt = spatialIndices.begin();
175 for ( auto it = mOverlayerSources.begin(); it != mOverlayerSources.end(); ++ it, ++spatialIteratorIt )
176 {
177 if ( feedback->isCanceled() )
178 break;
179
180 const QgsSpatialIndex &index = *spatialIteratorIt;
181 const QList<QgsFeatureId> matches = index.intersects( inputGeom.boundingBox() );
182 QVector< QgsGeometry > intersectingGeoms;
183 intersectingGeoms.reserve( matches.count() );
184 for ( const QgsFeatureId match : matches )
185 {
186 if ( feedback->isCanceled() )
187 break;
188
189 const QgsGeometry overlayGeometry = index.geometry( match );
190 if ( bufferGeomEngine->intersects( overlayGeometry.constGet() ) )
191 {
192 intersectingGeoms.append( overlayGeometry );
193 }
194 }
195
196 if ( feedback->isCanceled() )
197 break;
198
199 // dissolve intersecting features, calculate total area of them within our buffer
200 const QgsGeometry overlayDissolved = QgsGeometry::unaryUnion( intersectingGeoms, geometryParameters );
201
202 if ( feedback->isCanceled() )
203 break;
204
205 const QgsGeometry overlayIntersection = inputGeom.intersection( overlayDissolved, geometryParameters );
206
207 double overlayArea = 0;
208 try
209 {
210 overlayArea = da.measureArea( overlayIntersection );
211 }
212 catch ( QgsCsException & )
213 {
214 throw QgsProcessingException( QObject::tr( "An error occurred while calculating feature area" ) );
215 }
216
217 outAttributes.append( overlayArea );
218 outAttributes.append( 100 * overlayArea / inputArea );
219 }
220 }
221 else
222 {
223 // input feature has no geometry
224 for ( auto it = mOverlayerSources.begin(); it != mOverlayerSources.end(); ++ it )
225 {
226 outAttributes.append( QVariant() );
227 outAttributes.append( QVariant() );
228 }
229 }
230
231 feature.setAttributes( outAttributes );
232 if ( !sink->addFeature( feature, QgsFeatureSink::FastInsert ) )
233 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
234
235 i++;
236 feedback->setProgress( i * step );
237 }
238
239 QVariantMap outputs;
240 outputs.insert( QStringLiteral( "OUTPUT" ), destId );
241 return outputs;
242}
243
@ 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.
Custom exception class for Coordinate Reference System related exceptions.
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 & setSubsetOfAttributes(const QgsAttributeList &attrs)
Set a subset of attributes that will be fetched.
QgsFeatureRequest & setDestinationCrs(const QgsCoordinateReferenceSystem &crs, const QgsCoordinateTransformContext &context)
Sets the destination crs for feature's geometries.
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:58
QgsAttributes attributes
Definition qgsfeature.h:67
void setAttributes(const QgsAttributes &attrs)
Sets the feature's attributes.
QgsGeometry geometry
Definition qgsfeature.h:69
bool hasGeometry() const
Returns true if the feature has an associated geometry.
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.
void setGridSize(double size)
Sets the grid size which will be used to snap vertices of a geometry.
A geometry is the spatial representation of a feature.
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:76
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.
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:5917
qint64 QgsFeatureId
64 bit feature ids negative numbers are used for uncommitted/newly added features
QList< int > QgsAttributeList
Definition qgsfield.h:27