QGIS API Documentation  3.2.0-Bonn (bc43194)
qgsoverlayutils.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsoverlayutils.cpp
3  ---------------------
4  Date : April 2018
5  Copyright : (C) 2018 by Martin Dobias
6  Email : wonder dot sk at gmail dot com
7  ***************************************************************************
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  ***************************************************************************/
15 
16 #include "qgsoverlayutils.h"
17 
18 #include "qgsgeometryengine.h"
19 #include "qgsprocessingalgorithm.h"
20 
22 
24 static bool sanitizeIntersectionResult( QgsGeometry &geom, QgsWkbTypes::GeometryType geometryType )
25 {
26  if ( geom.isNull() )
27  {
28  // TODO: not sure if this ever happens - if it does, that means GEOS failed badly - would be good to have a test for such situation
29  throw QgsProcessingException( QStringLiteral( "%1\n\n%2" ).arg( QObject::tr( "GEOS geoprocessing error: intersection failed." ), geom.lastError() ) );
30  }
31 
32  // Intersection of geometries may give use also geometries we do not want in our results.
33  // For example, two square polygons touching at the corner have a point as the intersection, but no area.
34  // In other cases we may get a mixture of geometries in the output - we want to keep only the expected types.
36  {
37  // try to filter out irrelevant parts with different geometry type than what we want
38  geom.convertGeometryCollectionToSubclass( geometryType );
39  if ( geom.isEmpty() )
40  return false;
41  }
42 
43  if ( QgsWkbTypes::geometryType( geom.wkbType() ) != geometryType )
44  {
45  // we can't make use of this resulting geometry
46  return false;
47  }
48 
49  // some data providers are picky about the geometries we pass to them: we can't add single-part geometries
50  // when we promised multi-part geometries, so ensure we have the right type
51  geom.convertToMultiType();
52 
53  return true;
54 }
55 
56 
58 static bool sanitizeDifferenceResult( QgsGeometry &geom )
59 {
60  if ( geom.isNull() )
61  {
62  // TODO: not sure if this ever happens - if it does, that means GEOS failed badly - would be good to have a test for such situation
63  throw QgsProcessingException( QStringLiteral( "%1\n\n%2" ).arg( QObject::tr( "GEOS geoprocessing error: difference failed." ), geom.lastError() ) );
64  }
65 
66  // if geomB covers the whole source geometry, we get an empty geometry collection
67  if ( geom.isEmpty() )
68  return false;
69 
70  // some data providers are picky about the geometries we pass to them: we can't add single-part geometries
71  // when we promised multi-part geometries, so ensure we have the right type
72  geom.convertToMultiType();
73 
74  return true;
75 }
76 
77 
78 void QgsOverlayUtils::difference( const QgsFeatureSource &sourceA, const QgsFeatureSource &sourceB, QgsFeatureSink &sink, QgsProcessingContext &context, QgsProcessingFeedback *feedback, int &count, int totalCount, QgsOverlayUtils::DifferenceOutput outputAttrs )
79 {
80  QgsFeatureRequest requestB;
82  if ( outputAttrs != OutputBA )
83  requestB.setDestinationCrs( sourceA.sourceCrs(), context.transformContext() );
84  QgsSpatialIndex indexB( sourceB.getFeatures( requestB ), feedback );
85 
86  int fieldsCountA = sourceA.fields().count();
87  int fieldsCountB = sourceB.fields().count();
88  QgsAttributes attrs;
89  attrs.resize( outputAttrs == OutputA ? fieldsCountA : ( fieldsCountA + fieldsCountB ) );
90 
91  if ( totalCount == 0 )
92  totalCount = 1; // avoid division by zero
93 
94  QgsFeature featA;
95  QgsFeatureRequest requestA;
96  if ( outputAttrs == OutputBA )
97  requestA.setDestinationCrs( sourceB.sourceCrs(), context.transformContext() );
98  QgsFeatureIterator fitA = sourceA.getFeatures( requestA );
99  while ( fitA.nextFeature( featA ) )
100  {
101  if ( feedback->isCanceled() )
102  break;
103 
104  if ( featA.hasGeometry() )
105  {
106  QgsGeometry geom( featA.geometry() );
107  QgsFeatureIds intersects = indexB.intersects( geom.boundingBox() ).toSet();
108 
109  QgsFeatureRequest request;
110  request.setFilterFids( intersects );
112  if ( outputAttrs != OutputBA )
113  request.setDestinationCrs( sourceA.sourceCrs(), context.transformContext() );
114 
115  std::unique_ptr< QgsGeometryEngine > engine;
116  if ( !intersects.isEmpty() )
117  {
118  // use prepared geometries for faster intersection tests
119  engine.reset( QgsGeometry::createGeometryEngine( geom.constGet() ) );
120  engine->prepareGeometry();
121  }
122 
123  QVector<QgsGeometry> geometriesB;
124  QgsFeature featB;
125  QgsFeatureIterator fitB = sourceB.getFeatures( request );
126  while ( fitB.nextFeature( featB ) )
127  {
128  if ( feedback->isCanceled() )
129  break;
130 
131  if ( engine->intersects( featB.geometry().constGet() ) )
132  geometriesB << featB.geometry();
133  }
134 
135  if ( !geometriesB.isEmpty() )
136  {
137  QgsGeometry geomB = QgsGeometry::unaryUnion( geometriesB );
138  geom = geom.difference( geomB );
139  }
140 
141  if ( !sanitizeDifferenceResult( geom ) )
142  continue;
143 
144  const QgsAttributes attrsA( featA.attributes() );
145  switch ( outputAttrs )
146  {
147  case OutputA:
148  attrs = attrsA;
149  break;
150  case OutputAB:
151  for ( int i = 0; i < fieldsCountA; ++i )
152  attrs[i] = attrsA[i];
153  break;
154  case OutputBA:
155  for ( int i = 0; i < fieldsCountA; ++i )
156  attrs[i + fieldsCountB] = attrsA[i];
157  break;
158  }
159 
160  QgsFeature outFeat;
161  outFeat.setGeometry( geom );
162  outFeat.setAttributes( attrs );
163  sink.addFeature( outFeat, QgsFeatureSink::FastInsert );
164  }
165  else
166  {
167  // TODO: should we write out features that do not have geometry?
168  sink.addFeature( featA, QgsFeatureSink::FastInsert );
169  }
170 
171  ++count;
172  feedback->setProgress( count / ( double ) totalCount * 100. );
173  }
174 }
175 
176 
177 void QgsOverlayUtils::intersection( const QgsFeatureSource &sourceA, const QgsFeatureSource &sourceB, QgsFeatureSink &sink, QgsProcessingContext &context, QgsProcessingFeedback *feedback, int &count, int totalCount, const QList<int> &fieldIndicesA, const QList<int> &fieldIndicesB )
178 {
180  int attrCount = fieldIndicesA.count() + fieldIndicesB.count();
181 
182  QgsFeatureRequest request;
184  request.setDestinationCrs( sourceA.sourceCrs(), context.transformContext() );
185 
186  QgsFeature outFeat;
187  QgsSpatialIndex indexB( sourceB.getFeatures( request ), feedback );
188 
189  if ( totalCount == 0 )
190  totalCount = 1; // avoid division by zero
191 
192  QgsFeature featA;
193  QgsFeatureIterator fitA = sourceA.getFeatures( QgsFeatureRequest().setSubsetOfAttributes( fieldIndicesA ) );
194  while ( fitA.nextFeature( featA ) )
195  {
196  if ( feedback->isCanceled() )
197  break;
198 
199  if ( !featA.hasGeometry() )
200  continue;
201 
202  QgsGeometry geom( featA.geometry() );
203  QgsFeatureIds intersects = indexB.intersects( geom.boundingBox() ).toSet();
204 
205  QgsFeatureRequest request;
206  request.setFilterFids( intersects );
207  request.setDestinationCrs( sourceA.sourceCrs(), context.transformContext() );
208  request.setSubsetOfAttributes( fieldIndicesB );
209 
210  std::unique_ptr< QgsGeometryEngine > engine;
211  if ( !intersects.isEmpty() )
212  {
213  // use prepared geometries for faster intersection tests
214  engine.reset( QgsGeometry::createGeometryEngine( geom.constGet() ) );
215  engine->prepareGeometry();
216  }
217 
218  QgsAttributes outAttributes( attrCount );
219  const QgsAttributes attrsA( featA.attributes() );
220  for ( int i = 0; i < fieldIndicesA.count(); ++i )
221  outAttributes[i] = attrsA[fieldIndicesA[i]];
222 
223  QgsFeature featB;
224  QgsFeatureIterator fitB = sourceB.getFeatures( request );
225  while ( fitB.nextFeature( featB ) )
226  {
227  if ( feedback->isCanceled() )
228  break;
229 
230  QgsGeometry tmpGeom( featB.geometry() );
231  if ( !engine->intersects( tmpGeom.constGet() ) )
232  continue;
233 
234  QgsGeometry intGeom = geom.intersection( tmpGeom );
235  if ( !sanitizeIntersectionResult( intGeom, geometryType ) )
236  continue;
237 
238  const QgsAttributes attrsB( featB.attributes() );
239  for ( int i = 0; i < fieldIndicesB.count(); ++i )
240  outAttributes[fieldIndicesA.count() + i] = attrsB[fieldIndicesB[i]];
241 
242  outFeat.setGeometry( intGeom );
243  outFeat.setAttributes( outAttributes );
244  sink.addFeature( outFeat, QgsFeatureSink::FastInsert );
245  }
246 
247  ++count;
248  feedback->setProgress( count / ( double ) totalCount * 100. );
249  }
250 }
251 
252 void QgsOverlayUtils::resolveOverlaps( const QgsFeatureSource &source, QgsFeatureSink &sink, QgsProcessingFeedback *feedback )
253 {
254  int count = 0;
255  int totalCount = source.featureCount();
256  if ( totalCount == 0 )
257  return; // nothing to do here
258 
259  QgsFeatureId newFid = -1;
260 
262 
263  QgsFeatureRequest requestOnlyGeoms;
264  requestOnlyGeoms.setSubsetOfAttributes( QgsAttributeList() );
265 
266  QgsFeatureRequest requestOnlyAttrs;
267  requestOnlyAttrs.setFlags( QgsFeatureRequest::NoGeometry );
268 
269  QgsFeatureRequest requestOnlyIds;
270  requestOnlyIds.setFlags( QgsFeatureRequest::NoGeometry );
271  requestOnlyIds.setSubsetOfAttributes( QgsAttributeList() );
272 
273  // make a set of used feature IDs so that we do not try to reuse them for newly added features
274  QgsFeature f;
275  QSet<QgsFeatureId> fids;
276  QgsFeatureIterator it = source.getFeatures( requestOnlyIds );
277  while ( it.nextFeature( f ) )
278  {
279  if ( feedback->isCanceled() )
280  return;
281 
282  fids.insert( f.id() );
283  }
284 
285  QHash<QgsFeatureId, QgsGeometry> geometries;
286  QgsSpatialIndex index;
287  QHash<QgsFeatureId, QList<QgsFeatureId> > intersectingIds; // which features overlap a particular area
288 
289  // resolve intersections
290 
291  it = source.getFeatures( requestOnlyGeoms );
292  while ( it.nextFeature( f ) )
293  {
294  if ( feedback->isCanceled() )
295  return;
296 
297  QgsFeatureId fid1 = f.id();
298  QgsGeometry g1 = f.geometry();
299  std::unique_ptr< QgsGeometryEngine > g1engine;
300 
301  geometries.insert( fid1, g1 );
302  index.insertFeature( f );
303 
304  QgsRectangle bbox( f.geometry().boundingBox() );
305  const QList<QgsFeatureId> ids = index.intersects( bbox );
306  for ( QgsFeatureId fid2 : ids )
307  {
308  if ( fid1 == fid2 )
309  continue;
310 
311  if ( !g1engine )
312  {
313  // use prepared geometries for faster intersection tests
314  g1engine.reset( QgsGeometry::createGeometryEngine( g1.constGet() ) );
315  g1engine->prepareGeometry();
316  }
317 
318  QgsGeometry g2 = geometries.value( fid2 );
319  if ( !g1engine->intersects( g2.constGet() ) )
320  continue;
321 
322  QgsGeometry geomIntersection = g1.intersection( g2 );
323  if ( !sanitizeIntersectionResult( geomIntersection, geometryType ) )
324  continue;
325 
326  //
327  // add intersection geometry
328  //
329 
330  // figure out new fid
331  while ( fids.contains( newFid ) )
332  --newFid;
333  fids.insert( newFid );
334 
335  geometries.insert( newFid, geomIntersection );
336  QgsFeature fx( newFid );
337  fx.setGeometry( geomIntersection );
338 
339  index.insertFeature( fx );
340 
341  // figure out which feature IDs belong to this intersection. Some of the IDs can be of the newly
342  // created geometries - in such case we need to retrieve original IDs
343  QList<QgsFeatureId> lst;
344  if ( intersectingIds.contains( fid1 ) )
345  lst << intersectingIds.value( fid1 );
346  else
347  lst << fid1;
348  if ( intersectingIds.contains( fid2 ) )
349  lst << intersectingIds.value( fid2 );
350  else
351  lst << fid2;
352  intersectingIds.insert( newFid, lst );
353 
354  //
355  // update f1
356  //
357 
358  QgsGeometry g12 = g1.difference( g2 );
359 
360  index.deleteFeature( f );
361  geometries.remove( fid1 );
362 
363  if ( sanitizeDifferenceResult( g12 ) )
364  {
365  geometries.insert( fid1, g12 );
366 
367  QgsFeature f1x( fid1 );
368  f1x.setGeometry( g12 );
369  index.insertFeature( f1x );
370  }
371 
372  //
373  // update f2
374  //
375 
376  QgsGeometry g21 = g2.difference( g1 );
377 
378  QgsFeature f2old( fid2 );
379  f2old.setGeometry( g2 );
380  index.deleteFeature( f2old );
381 
382  geometries.remove( fid2 );
383 
384  if ( sanitizeDifferenceResult( g21 ) )
385  {
386  geometries.insert( fid2, g21 );
387 
388  QgsFeature f2x( fid2 );
389  f2x.setGeometry( g21 );
390  index.insertFeature( f2x );
391  }
392 
393  // update our temporary copy of the geometry to what is left from it
394  g1 = g12;
395  g1engine.reset();
396  }
397 
398  ++count;
399  feedback->setProgress( count / ( double ) totalCount * 100. );
400  }
401 
402  // release some memory of structures we don't need anymore
403 
404  fids.clear();
405  index = QgsSpatialIndex();
406 
407  // load attributes
408 
409  QHash<QgsFeatureId, QgsAttributes> attributesHash;
410  it = source.getFeatures( requestOnlyAttrs );
411  while ( it.nextFeature( f ) )
412  {
413  if ( feedback->isCanceled() )
414  return;
415 
416  attributesHash.insert( f.id(), f.attributes() );
417  }
418 
419  // store stuff in the sink
420 
421  for ( auto i = geometries.constBegin(); i != geometries.constEnd(); ++i )
422  {
423  if ( feedback->isCanceled() )
424  return;
425 
426  QgsFeature outFeature( i.key() );
427  outFeature.setGeometry( i.value() );
428 
429  if ( intersectingIds.contains( i.key() ) )
430  {
431  const QList<QgsFeatureId> ids = intersectingIds.value( i.key() );
432  for ( QgsFeatureId id : ids )
433  {
434  outFeature.setAttributes( attributesHash.value( id ) );
435  sink.addFeature( outFeature, QgsFeatureSink::FastInsert );
436  }
437  }
438  else
439  {
440  outFeature.setAttributes( attributesHash.value( i.key() ) );
441  sink.addFeature( outFeature, QgsFeatureSink::FastInsert );
442  }
443  }
444 }
445 
bool convertGeometryCollectionToSubclass(QgsWkbTypes::GeometryType geomType)
Converts geometry collection to a the desired geometry type subclass (multi-point, multi-linestring or multi-polygon).
QgsFeatureRequest & setDestinationCrs(const QgsCoordinateReferenceSystem &crs, const QgsCoordinateTransformContext &context)
Sets the destination crs for feature&#39;s geometries.
QgsFeatureId id
Definition: qgsfeature.h:71
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:40
static Type multiType(Type type)
Returns the multi type for a WKB type.
Definition: qgswkbtypes.h:296
Base class for providing feedback from a processing algorithm.
virtual QgsFields fields() const =0
Returns the fields associated with features in the source.
bool isNull() const
Returns true if the geometry is null (ie, contains no underlying geometry accessible via geometry() )...
QSet< QgsFeatureId > QgsFeatureIds
Definition: qgsfeature.h:544
QgsWkbTypes::Type wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
virtual QgsWkbTypes::Type wkbType() const =0
Returns the geometry type for features returned by this source.
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
An interface for objects which accept features via addFeature(s) methods.
QgsFeatureRequest & setSubsetOfAttributes(const QgsAttributeList &attrs)
Set a subset of attributes that will be fetched.
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:104
void setAttributes(const QgsAttributes &attrs)
Sets the feature&#39;s attributes.
Definition: qgsfeature.cpp:127
bool insertFeature(const QgsFeature &feature)
Adds a feature to the index.
QgsGeometry intersection(const QgsGeometry &geometry) const
Returns a geometry representing the points shared by this geometry and other.
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
Definition: qgsfeature.h:62
bool hasGeometry() const
Returns true if the feature has an associated geometry.
Definition: qgsfeature.cpp:190
int count() const
Returns number of items.
Definition: qgsfields.cpp:115
bool intersects(const QgsRectangle &rectangle) const
Returns true if this geometry exactly intersects with a rectangle.
static GeometryType geometryType(Type type)
Returns the geometry type for a WKB type, e.g., both MultiPolygon and CurvePolygon would have a Polyg...
Definition: qgswkbtypes.h:663
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 isEmpty() const
Returns true if the geometry is empty (eg a linestring with no vertices, or a collection with no geom...
virtual bool addFeature(QgsFeature &feature, QgsFeatureSink::Flags flags=nullptr)
Adds a single feature to the sink.
QgsGeometry geometry() const
Returns the geometry associated with this feature.
Definition: qgsfeature.cpp:101
virtual QgsCoordinateReferenceSystem sourceCrs() const =0
Returns the coordinate reference system for features in the source.
static QgsGeometryEngine * createGeometryEngine(const QgsAbstractGeometry *geometry)
Creates and returns a new geometry engine.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
GeometryType
The geometry types are used to group QgsWkbTypes::Type in a coarse way.
Definition: qgswkbtypes.h:137
bool deleteFeature(const QgsFeature &feature)
Removes a feature from the index.
QString lastError() const
Returns an error string referring to the last error encountered either when this geometry was created...
A spatial index for QgsFeature objects.
bool convertToMultiType()
Converts single type geometry into multitype geometry e.g.
QgsFeatureRequest & setFilterFids(const QgsFeatureIds &fids)
Sets feature IDs that should be fetched.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
An interface for objects which provide features via a getFeatures method.
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
qint64 QgsFeatureId
Definition: qgsfeature.h:37
QList< int > QgsAttributeList
Definition: qgsfield.h:27
bool nextFeature(QgsFeature &f)
Geometry is not required. It may still be returned if e.g. required for a filter condition.
A vector of attributes.
Definition: qgsattributes.h:58
static Type flatType(Type type)
Returns the flat type for a WKB type.
Definition: qgswkbtypes.h:427
Contains information about the context in which a processing algorithm is executed.
static QgsGeometry unaryUnion(const QVector< QgsGeometry > &geometries)
Compute the unary union on a list of geometries.
QList< QgsFeatureId > intersects(const QgsRectangle &rectangle) const
Returns a list of features with a bounding box which intersects the specified rectangle.
QgsAttributes attributes
Definition: qgsfeature.h:72
virtual QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest()) const =0
Returns an iterator for the features in the source.
QgsFeatureRequest & setFlags(QgsFeatureRequest::Flags flags)
Sets flags that affect how features will be fetched.
virtual long featureCount() const =0
Returns the number of features contained in the source, or -1 if the feature count is unknown...
QgsGeometry difference(const QgsGeometry &geometry) const
Returns a geometry representing the points making up this geometry that do not make up other...