QGIS API Documentation  3.14.0-Pi (9f7028fd23)
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 
23 
24 QString QgsCalculateVectorOverlapsAlgorithm::name() const
25 {
26  return QStringLiteral( "calculatevectoroverlaps" );
27 }
28 
29 QString QgsCalculateVectorOverlapsAlgorithm::displayName() const
30 {
31  return QObject::tr( "Overlap analysis" );
32 }
33 
34 QStringList QgsCalculateVectorOverlapsAlgorithm::tags() const
35 {
36  return QObject::tr( "vector,overlay,area,percentage,intersection" ).split( ',' );
37 }
38 
39 QString QgsCalculateVectorOverlapsAlgorithm::group() const
40 {
41  return QObject::tr( "Vector analysis" );
42 }
43 
44 QString QgsCalculateVectorOverlapsAlgorithm::groupId() const
45 {
46  return QStringLiteral( "vectoranalysis" );
47 }
48 
49 void QgsCalculateVectorOverlapsAlgorithm::initAlgorithm( const QVariantMap & )
50 {
51  addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT" ), QObject::tr( "Input layer" ), QList< int >() << QgsProcessing::TypeVectorPolygon ) );
52  addParameter( new QgsProcessingParameterMultipleLayers( QStringLiteral( "LAYERS" ), QObject::tr( "Overlay layers" ), QgsProcessing::TypeVectorPolygon ) );
53  addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Overlap" ) ) );
54 }
55 
56 QIcon QgsCalculateVectorOverlapsAlgorithm::icon() const
57 {
58  return QgsApplication::getThemeIcon( QStringLiteral( "/algorithms/mAlgorithmClip.svg" ) );
59 }
60 
61 QString QgsCalculateVectorOverlapsAlgorithm::svgIconPath() const
62 {
63  return QgsApplication::iconPath( QStringLiteral( "/algorithms/mAlgorithmClip.svg" ) );
64 }
65 
66 QString QgsCalculateVectorOverlapsAlgorithm::shortHelpString() const
67 {
68  return QObject::tr( "This algorithm calculates the area and percentage cover by which features from an input layer "
69  "are overlapped by features from a selection of overlay layers.\n\n"
70  "New attributes are added to the output layer reporting the total area of overlap and percentage of the input feature overlapped "
71  "by each of the selected overlay layers." );
72 }
73 
74 QgsCalculateVectorOverlapsAlgorithm *QgsCalculateVectorOverlapsAlgorithm::createInstance() const
75 {
76  return new QgsCalculateVectorOverlapsAlgorithm();
77 }
78 
79 bool QgsCalculateVectorOverlapsAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
80 {
81  mSource.reset( parameterAsSource( parameters, QStringLiteral( "INPUT" ), context ) );
82  if ( !mSource )
83  throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT" ) ) );
84 
85  mOutputFields = mSource->fields();
86 
87  const QList< QgsMapLayer * > layers = parameterAsLayerList( parameters, QStringLiteral( "LAYERS" ), context );
88  mOverlayerSources.reserve( layers.size() );
89  mLayerNames.reserve( layers.size() );
90  for ( QgsMapLayer *layer : layers )
91  {
92  if ( QgsVectorLayer *vl = qobject_cast< QgsVectorLayer * >( layer ) )
93  {
94  mLayerNames << layer->name();
95  mOverlayerSources.emplace_back( qgis::make_unique< QgsVectorLayerFeatureSource >( vl ) );
96  mOutputFields.append( QgsField( QStringLiteral( "%1_area" ).arg( vl->name() ), QVariant::Double ) );
97  mOutputFields.append( QgsField( QStringLiteral( "%1_pc" ).arg( vl->name() ), QVariant::Double ) );
98  }
99  }
100 
101  mOutputType = mSource->wkbType();
102  mCrs = mSource->sourceCrs();
103  mInputCount = mSource->featureCount();
104  mInputFeatures = mSource->getFeatures();
105  return true;
106 }
107 
108 QVariantMap QgsCalculateVectorOverlapsAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
109 {
110  QString destId;
111  std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, destId, mOutputFields,
112  mOutputType, mCrs ) );
113  if ( !sink )
114  throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );
115 
116  // build a spatial index for each constraint layer for speed. We also store input constraint geometries here,
117  // to avoid refetching and projecting them later
118  QList< QgsSpatialIndex > spatialIndices;
119  spatialIndices.reserve( mLayerNames.size() );
120  auto nameIt = mLayerNames.constBegin();
121  for ( auto sourceIt = mOverlayerSources.begin(); sourceIt != mOverlayerSources.end(); ++sourceIt, ++nameIt )
122  {
123  feedback->pushInfo( QObject::tr( "Preparing %1" ).arg( *nameIt ) );
124  QgsFeatureIterator featureIt = ( *sourceIt )->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( QgsAttributeList() ).setDestinationCrs( mCrs, context.transformContext() ).setInvalidGeometryCheck( context.invalidGeometryCheck() ).setInvalidGeometryCallback( context.invalidGeometryCallback() ) );
125  spatialIndices << QgsSpatialIndex( featureIt, feedback, QgsSpatialIndex::FlagStoreFeatureGeometries );
126  }
127 
128  QgsDistanceArea da;
129  da.setSourceCrs( mCrs, context.transformContext() );
130  da.setEllipsoid( context.project()->ellipsoid() );
131 
132  // loop through input
133  double step = mInputCount > 0 ? 100.0 / mInputCount : 0;
134  long i = 0;
135  QgsFeature feature;
136  while ( mInputFeatures.nextFeature( feature ) )
137  {
138  if ( feedback->isCanceled() )
139  break;
140 
141  QgsAttributes outAttributes = feature.attributes();
142  if ( feature.hasGeometry() && !qgsDoubleNear( feature.geometry().area(), 0.0 ) )
143  {
144  const QgsGeometry inputGeom = feature.geometry();
145  const double inputArea = da.measureArea( inputGeom );
146 
147  // prepare for lots of intersection tests (for speed)
148  std::unique_ptr< QgsGeometryEngine > bufferGeomEngine( QgsGeometry::createGeometryEngine( inputGeom.constGet() ) );
149  bufferGeomEngine->prepareGeometry();
150 
151  // calculate overlap attributes
152  auto spatialIteratorIt = spatialIndices.begin();
153  for ( auto it = mOverlayerSources.begin(); it != mOverlayerSources.end(); ++ it, ++spatialIteratorIt )
154  {
155  if ( feedback->isCanceled() )
156  break;
157 
158  QgsSpatialIndex &index = *spatialIteratorIt;
159  const QList<QgsFeatureId> matches = index.intersects( inputGeom.boundingBox() );
160  QVector< QgsGeometry > intersectingGeoms;
161  intersectingGeoms.reserve( matches.count() );
162  for ( QgsFeatureId match : matches )
163  {
164  if ( feedback->isCanceled() )
165  break;
166 
167  const QgsGeometry overlayGeometry = index.geometry( match );
168  if ( bufferGeomEngine->intersects( overlayGeometry.constGet() ) )
169  {
170  intersectingGeoms.append( overlayGeometry );
171  }
172  }
173 
174  if ( feedback->isCanceled() )
175  break;
176 
177  // dissolve intersecting features, calculate total area of them within our buffer
178  QgsGeometry overlayDissolved = QgsGeometry::unaryUnion( intersectingGeoms );
179 
180  if ( feedback->isCanceled() )
181  break;
182 
183  QgsGeometry overlayIntersection = inputGeom.intersection( overlayDissolved );
184 
185  const double overlayArea = da.measureArea( overlayIntersection );
186  outAttributes.append( overlayArea );
187  outAttributes.append( 100 * overlayArea / inputArea );
188  }
189  }
190  else
191  {
192  // input feature has no geometry
193  for ( auto it = mOverlayerSources.begin(); it != mOverlayerSources.end(); ++ it )
194  {
195  outAttributes.append( QVariant() );
196  outAttributes.append( QVariant() );
197  }
198  }
199 
200  feature.setAttributes( outAttributes );
201  sink->addFeature( feature, QgsFeatureSink::FastInsert );
202 
203  i++;
204  feedback->setProgress( i * step );
205  }
206 
207  QVariantMap outputs;
208  outputs.insert( QStringLiteral( "OUTPUT" ), destId );
209  return outputs;
210 }
211 
QgsFeedback::setProgress
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:75
QgsFeatureRequest::setInvalidGeometryCheck
QgsFeatureRequest & setInvalidGeometryCheck(InvalidGeometryCheck check)
Sets invalid geometry checking behavior.
Definition: qgsfeaturerequest.cpp:117
QgsApplication::getThemeIcon
static QIcon getThemeIcon(const QString &name)
Helper to get a theme icon.
Definition: qgsapplication.cpp:605
QgsProcessingContext::project
QgsProject * project() const
Returns the project in which the algorithm is being executed.
Definition: qgsprocessingcontext.h:99
QgsProcessingFeedback
Definition: qgsprocessingfeedback.h:37
QgsProcessingFeedback::pushInfo
virtual void pushInfo(const QString &info)
Pushes a general informational message from the algorithm.
Definition: qgsprocessingfeedback.cpp:48
qgsalgorithmcalculateoverlaps.h
QgsProcessing::TypeVectorPolygon
@ TypeVectorPolygon
Vector polygon layers.
Definition: qgsprocessing.h:50
QgsFeature::geometry
QgsGeometry geometry
Definition: qgsfeature.h:71
QgsDistanceArea::measureArea
double measureArea(const QgsGeometry &geometry) const
Measures the area of a geometry.
Definition: qgsdistancearea.cpp:183
QgsFeatureRequest::setInvalidGeometryCallback
QgsFeatureRequest & setInvalidGeometryCallback(const std::function< void(const QgsFeature &)> &callback)
Sets a callback function to use when encountering an invalid geometry and invalidGeometryCheck() is s...
Definition: qgsfeaturerequest.cpp:123
QgsProcessingParameterFeatureSource
Definition: qgsprocessingparameters.h:2612
QgsApplication::iconPath
static QString iconPath(const QString &iconFile)
Returns path to the desired icon file.
Definition: qgsapplication.cpp:594
QgsSpatialIndex::geometry
QgsGeometry geometry(QgsFeatureId id) const
Returns the stored geometry for the indexed feature with matching id.
Definition: qgsspatialindex.cpp:552
QgsAttributeList
QList< int > QgsAttributeList
Definition: qgsfield.h:26
QgsProcessingParameterMultipleLayers
Definition: qgsprocessingparameters.h:1756
QgsGeometry::intersection
QgsGeometry intersection(const QgsGeometry &geometry) const
Returns a geometry representing the points shared by this geometry and other.
Definition: qgsgeometry.cpp:2395
QgsProcessingParameterFeatureSink
Definition: qgsprocessingparameters.h:2773
QgsDistanceArea::setEllipsoid
bool setEllipsoid(const QString &ellipsoid)
Sets the ellipsoid by its acronym.
Definition: qgsdistancearea.cpp:66
QgsFeatureRequest
Definition: qgsfeaturerequest.h:75
qgsgeometryengine.h
QgsProcessingContext::invalidGeometryCheck
QgsFeatureRequest::InvalidGeometryCheck invalidGeometryCheck() const
Returns the behavior used for checking invalid geometries in input layers.
Definition: qgsprocessingcontext.h:296
QgsProcessingContext
Definition: qgsprocessingcontext.h:43
qgsDoubleNear
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:315
QgsSpatialIndex
Definition: qgsspatialindex.h:67
QgsProcessingContext::transformContext
QgsCoordinateTransformContext transformContext() const
Returns the coordinate transform context.
Definition: qgsprocessingcontext.h:135
QgsGeometry::constGet
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
Definition: qgsgeometry.cpp:128
QgsDistanceArea::setSourceCrs
void setSourceCrs(const QgsCoordinateReferenceSystem &crs, const QgsCoordinateTransformContext &context)
Sets source spatial reference system crs.
Definition: qgsdistancearea.cpp:60
QgsFeature::attributes
QgsAttributes attributes
Definition: qgsfeature.h:69
qgsvectorlayer.h
QgsSpatialIndex::FlagStoreFeatureGeometries
@ FlagStoreFeatureGeometries
Indicates that the spatial index should also store feature geometries. This requires more memory,...
Definition: qgsspatialindex.h:77
QgsFeedback::isCanceled
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:66
QgsGeometry::createGeometryEngine
static QgsGeometryEngine * createGeometryEngine(const QgsAbstractGeometry *geometry)
Creates and returns a new geometry engine.
Definition: qgsgeometry.cpp:3659
QgsGeometry
Definition: qgsgeometry.h:122
QgsVectorLayer
Definition: qgsvectorlayer.h:385
QgsFeature::hasGeometry
bool hasGeometry() const
Returns true if the feature has an associated geometry.
Definition: qgsfeature.cpp:197
QgsMapLayer
Definition: qgsmaplayer.h:81
QgsSpatialIndex::intersects
QList< QgsFeatureId > intersects(const QgsRectangle &rectangle) const
Returns a list of features with a bounding box which intersects the specified rectangle.
Definition: qgsspatialindex.cpp:488
QgsGeometry::boundingBox
QgsRectangle boundingBox() const
Returns the bounding box of the geometry.
Definition: qgsgeometry.cpp:962
QgsDistanceArea
Definition: qgsdistancearea.h:49
QgsAttributes
Definition: qgsattributes.h:57
QgsFeature
Definition: qgsfeature.h:55
QgsFeature::setAttributes
void setAttributes(const QgsAttributes &attrs)
Sets the feature's attributes.
Definition: qgsfeature.cpp:127
QgsProcessingContext::invalidGeometryCallback
std::function< void(const QgsFeature &) > invalidGeometryCallback() const
Returns the callback function to use when encountering an invalid geometry and invalidGeometryCheck()...
Definition: qgsprocessingcontext.h:338
QgsFeatureIterator
Definition: qgsfeatureiterator.h:263
QgsGeometry::area
double area() const
Returns the planar, 2-dimensional area of the geometry.
Definition: qgsgeometry.cpp:1775
QgsProject::ellipsoid
QString ellipsoid
Definition: qgsproject.h:100
QgsProcessingException
Definition: qgsexception.h:82
QgsGeometry::unaryUnion
static QgsGeometry unaryUnion(const QVector< QgsGeometry > &geometries)
Compute the unary union on a list of geometries.
Definition: qgsgeometry.cpp:2785
QgsFeatureSink::FastInsert
@ FastInsert
Use faster inserts, at the cost of updating the passed features to reflect changes made at the provid...
Definition: qgsfeaturesink.h:70
QgsFeatureId
qint64 QgsFeatureId
Definition: qgsfeatureid.h:25
QgsField
Definition: qgsfield.h:49