QGIS API Documentation  3.17.0-Master (df2c9ff931)
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 
23 bool QgsOverlayUtils::sanitizeIntersectionResult( QgsGeometry &geom, QgsWkbTypes::GeometryType geometryType )
24 {
25  if ( geom.isNull() )
26  {
27  // 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
28  throw QgsProcessingException( QStringLiteral( "%1\n\n%2" ).arg( QObject::tr( "GEOS geoprocessing error: intersection failed." ), geom.lastError() ) );
29  }
30 
31  // Intersection of geometries may give use also geometries we do not want in our results.
32  // For example, two square polygons touching at the corner have a point as the intersection, but no area.
33  // In other cases we may get a mixture of geometries in the output - we want to keep only the expected types.
35  {
36  // try to filter out irrelevant parts with different geometry type than what we want
37  geom.convertGeometryCollectionToSubclass( geometryType );
38  if ( geom.isEmpty() )
39  return false;
40  }
41 
42  if ( QgsWkbTypes::geometryType( geom.wkbType() ) != geometryType )
43  {
44  // we can't make use of this resulting geometry
45  return false;
46  }
47 
48  // some data providers are picky about the geometries we pass to them: we can't add single-part geometries
49  // when we promised multi-part geometries, so ensure we have the right type
50  geom.convertToMultiType();
51 
52  return true;
53 }
54 
55 
57 static bool sanitizeDifferenceResult( QgsGeometry &geom, QgsWkbTypes::GeometryType geometryType )
58 {
59  if ( geom.isNull() )
60  {
61  // 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
62  throw QgsProcessingException( QStringLiteral( "%1\n\n%2" ).arg( QObject::tr( "GEOS geoprocessing error: difference failed." ), geom.lastError() ) );
63  }
64 
65  //fix geometry collections
67  {
68  // try to filter out irrelevant parts with different geometry type than what we want
69  geom.convertGeometryCollectionToSubclass( geometryType );
70  }
71 
72 
73  // if geomB covers the whole source geometry, we get an empty geometry collection
74  if ( geom.isEmpty() )
75  return false;
76 
77  // some data providers are picky about the geometries we pass to them: we can't add single-part geometries
78  // when we promised multi-part geometries, so ensure we have the right type
79  geom.convertToMultiType();
80 
81  return true;
82 }
83 
84 
85 void QgsOverlayUtils::difference( const QgsFeatureSource &sourceA, const QgsFeatureSource &sourceB, QgsFeatureSink &sink, QgsProcessingContext &context, QgsProcessingFeedback *feedback, int &count, int totalCount, QgsOverlayUtils::DifferenceOutput outputAttrs )
86 {
88  QgsFeatureRequest requestB;
89  requestB.setNoAttributes();
90  if ( outputAttrs != OutputBA )
91  requestB.setDestinationCrs( sourceA.sourceCrs(), context.transformContext() );
92  QgsSpatialIndex indexB( sourceB.getFeatures( requestB ), feedback );
93 
94  int fieldsCountA = sourceA.fields().count();
95  int fieldsCountB = sourceB.fields().count();
96  QgsAttributes attrs;
97  attrs.resize( outputAttrs == OutputA ? fieldsCountA : ( fieldsCountA + fieldsCountB ) );
98 
99  if ( totalCount == 0 )
100  totalCount = 1; // avoid division by zero
101 
102  QgsFeature featA;
103  QgsFeatureRequest requestA;
104  requestA.setInvalidGeometryCheck( context.invalidGeometryCheck() );
105  if ( outputAttrs == OutputBA )
106  requestA.setDestinationCrs( sourceB.sourceCrs(), context.transformContext() );
107  QgsFeatureIterator fitA = sourceA.getFeatures( requestA );
108  while ( fitA.nextFeature( featA ) )
109  {
110  if ( feedback->isCanceled() )
111  break;
112 
113  if ( featA.hasGeometry() )
114  {
115  QgsGeometry geom( featA.geometry() );
116  QgsFeatureIds intersects = qgis::listToSet( indexB.intersects( geom.boundingBox() ) );
117 
118  QgsFeatureRequest request;
119  request.setFilterFids( intersects );
120  request.setNoAttributes();
121  if ( outputAttrs != OutputBA )
122  request.setDestinationCrs( sourceA.sourceCrs(), context.transformContext() );
123 
124  std::unique_ptr< QgsGeometryEngine > engine;
125  if ( !intersects.isEmpty() )
126  {
127  // use prepared geometries for faster intersection tests
128  engine.reset( QgsGeometry::createGeometryEngine( geom.constGet() ) );
129  engine->prepareGeometry();
130  }
131 
132  QVector<QgsGeometry> geometriesB;
133  QgsFeature featB;
134  QgsFeatureIterator fitB = sourceB.getFeatures( request );
135  while ( fitB.nextFeature( featB ) )
136  {
137  if ( feedback->isCanceled() )
138  break;
139 
140  if ( engine->intersects( featB.geometry().constGet() ) )
141  geometriesB << featB.geometry();
142  }
143 
144  if ( !geometriesB.isEmpty() )
145  {
146  QgsGeometry geomB = QgsGeometry::unaryUnion( geometriesB );
147  if ( !geomB.lastError().isEmpty() )
148  {
149  // This may happen if input geometries from a layer do not line up well (for example polygons
150  // that are nearly touching each other, but there is a very tiny overlap or gap at one of the edges).
151  // It is possible to get rid of this issue in two steps:
152  // 1. snap geometries with a small tolerance (e.g. 1cm) using QgsGeometrySnapperSingleSource
153  // 2. fix geometries (removes polygons collapsed to lines etc.) using MakeValid
154  throw QgsProcessingException( QStringLiteral( "%1\n\n%2" ).arg( QObject::tr( "GEOS geoprocessing error: unary union failed." ), geomB.lastError() ) );
155  }
156  geom = geom.difference( geomB );
157  }
158 
159  if ( !sanitizeDifferenceResult( geom, geometryType ) )
160  continue;
161 
162  const QgsAttributes attrsA( featA.attributes() );
163  switch ( outputAttrs )
164  {
165  case OutputA:
166  attrs = attrsA;
167  break;
168  case OutputAB:
169  for ( int i = 0; i < fieldsCountA; ++i )
170  attrs[i] = attrsA[i];
171  break;
172  case OutputBA:
173  for ( int i = 0; i < fieldsCountA; ++i )
174  attrs[i + fieldsCountB] = attrsA[i];
175  break;
176  }
177 
178  QgsFeature outFeat;
179  outFeat.setGeometry( geom );
180  outFeat.setAttributes( attrs );
181  sink.addFeature( outFeat, QgsFeatureSink::FastInsert );
182  }
183  else
184  {
185  // TODO: should we write out features that do not have geometry?
186  sink.addFeature( featA, QgsFeatureSink::FastInsert );
187  }
188 
189  ++count;
190  feedback->setProgress( count / ( double ) totalCount * 100. );
191  }
192 }
193 
194 
195 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 )
196 {
198  int attrCount = fieldIndicesA.count() + fieldIndicesB.count();
199 
200  QgsFeatureRequest request;
201  request.setNoAttributes();
202  request.setDestinationCrs( sourceA.sourceCrs(), context.transformContext() );
203 
204  QgsFeature outFeat;
205  QgsSpatialIndex indexB( sourceB.getFeatures( request ), feedback );
206 
207  if ( totalCount == 0 )
208  totalCount = 1; // avoid division by zero
209 
210  QgsFeature featA;
211  QgsFeatureIterator fitA = sourceA.getFeatures( QgsFeatureRequest().setSubsetOfAttributes( fieldIndicesA ) );
212  while ( fitA.nextFeature( featA ) )
213  {
214  if ( feedback->isCanceled() )
215  break;
216 
217  if ( !featA.hasGeometry() )
218  continue;
219 
220  QgsGeometry geom( featA.geometry() );
221  QgsFeatureIds intersects = qgis::listToSet( indexB.intersects( geom.boundingBox() ) );
222 
223  QgsFeatureRequest request;
224  request.setFilterFids( intersects );
225  request.setDestinationCrs( sourceA.sourceCrs(), context.transformContext() );
226  request.setSubsetOfAttributes( fieldIndicesB );
227 
228  std::unique_ptr< QgsGeometryEngine > engine;
229  if ( !intersects.isEmpty() )
230  {
231  // use prepared geometries for faster intersection tests
232  engine.reset( QgsGeometry::createGeometryEngine( geom.constGet() ) );
233  engine->prepareGeometry();
234  }
235 
236  QgsAttributes outAttributes( attrCount );
237  const QgsAttributes attrsA( featA.attributes() );
238  for ( int i = 0; i < fieldIndicesA.count(); ++i )
239  outAttributes[i] = attrsA[fieldIndicesA[i]];
240 
241  QgsFeature featB;
242  QgsFeatureIterator fitB = sourceB.getFeatures( request );
243  while ( fitB.nextFeature( featB ) )
244  {
245  if ( feedback->isCanceled() )
246  break;
247 
248  QgsGeometry tmpGeom( featB.geometry() );
249  if ( !engine->intersects( tmpGeom.constGet() ) )
250  continue;
251 
252  QgsGeometry intGeom = geom.intersection( tmpGeom );
253  if ( !sanitizeIntersectionResult( intGeom, geometryType ) )
254  continue;
255 
256  const QgsAttributes attrsB( featB.attributes() );
257  for ( int i = 0; i < fieldIndicesB.count(); ++i )
258  outAttributes[fieldIndicesA.count() + i] = attrsB[fieldIndicesB[i]];
259 
260  outFeat.setGeometry( intGeom );
261  outFeat.setAttributes( outAttributes );
262  sink.addFeature( outFeat, QgsFeatureSink::FastInsert );
263  }
264 
265  ++count;
266  feedback->setProgress( count / ( double ) totalCount * 100. );
267  }
268 }
269 
270 void QgsOverlayUtils::resolveOverlaps( const QgsFeatureSource &source, QgsFeatureSink &sink, QgsProcessingFeedback *feedback )
271 {
272  int count = 0;
273  int totalCount = source.featureCount();
274  if ( totalCount == 0 )
275  return; // nothing to do here
276 
277  QgsFeatureId newFid = -1;
278 
280 
281  QgsFeatureRequest requestOnlyGeoms;
282  requestOnlyGeoms.setNoAttributes();
283 
284  QgsFeatureRequest requestOnlyAttrs;
285  requestOnlyAttrs.setFlags( QgsFeatureRequest::NoGeometry );
286 
287  QgsFeatureRequest requestOnlyIds;
288  requestOnlyIds.setFlags( QgsFeatureRequest::NoGeometry );
289  requestOnlyIds.setNoAttributes();
290 
291  // make a set of used feature IDs so that we do not try to reuse them for newly added features
292  QgsFeature f;
293  QSet<QgsFeatureId> fids;
294  QgsFeatureIterator it = source.getFeatures( requestOnlyIds );
295  while ( it.nextFeature( f ) )
296  {
297  if ( feedback->isCanceled() )
298  return;
299 
300  fids.insert( f.id() );
301  }
302 
303  QHash<QgsFeatureId, QgsGeometry> geometries;
304  QgsSpatialIndex index;
305  QHash<QgsFeatureId, QList<QgsFeatureId> > intersectingIds; // which features overlap a particular area
306 
307  // resolve intersections
308 
309  it = source.getFeatures( requestOnlyGeoms );
310  while ( it.nextFeature( f ) )
311  {
312  if ( feedback->isCanceled() )
313  return;
314 
315  QgsFeatureId fid1 = f.id();
316  QgsGeometry g1 = f.geometry();
317  std::unique_ptr< QgsGeometryEngine > g1engine;
318 
319  geometries.insert( fid1, g1 );
320  index.addFeature( f );
321 
322  QgsRectangle bbox( f.geometry().boundingBox() );
323  const QList<QgsFeatureId> ids = index.intersects( bbox );
324  for ( QgsFeatureId fid2 : ids )
325  {
326  if ( fid1 == fid2 )
327  continue;
328 
329  if ( !g1engine )
330  {
331  // use prepared geometries for faster intersection tests
332  g1engine.reset( QgsGeometry::createGeometryEngine( g1.constGet() ) );
333  g1engine->prepareGeometry();
334  }
335 
336  QgsGeometry g2 = geometries.value( fid2 );
337  if ( !g1engine->intersects( g2.constGet() ) )
338  continue;
339 
340  QgsGeometry geomIntersection = g1.intersection( g2 );
341  if ( !sanitizeIntersectionResult( geomIntersection, geometryType ) )
342  continue;
343 
344  //
345  // add intersection geometry
346  //
347 
348  // figure out new fid
349  while ( fids.contains( newFid ) )
350  --newFid;
351  fids.insert( newFid );
352 
353  geometries.insert( newFid, geomIntersection );
354  QgsFeature fx( newFid );
355  fx.setGeometry( geomIntersection );
356 
357  index.addFeature( fx );
358 
359  // figure out which feature IDs belong to this intersection. Some of the IDs can be of the newly
360  // created geometries - in such case we need to retrieve original IDs
361  QList<QgsFeatureId> lst;
362  if ( intersectingIds.contains( fid1 ) )
363  lst << intersectingIds.value( fid1 );
364  else
365  lst << fid1;
366  if ( intersectingIds.contains( fid2 ) )
367  lst << intersectingIds.value( fid2 );
368  else
369  lst << fid2;
370  intersectingIds.insert( newFid, lst );
371 
372  //
373  // update f1
374  //
375 
376  QgsGeometry g12 = g1.difference( g2 );
377 
378  index.deleteFeature( f );
379  geometries.remove( fid1 );
380 
381  if ( sanitizeDifferenceResult( g12, geometryType ) )
382  {
383  geometries.insert( fid1, g12 );
384 
385  QgsFeature f1x( fid1 );
386  f1x.setGeometry( g12 );
387  index.addFeature( f1x );
388  }
389 
390  //
391  // update f2
392  //
393 
394  QgsGeometry g21 = g2.difference( g1 );
395 
396  QgsFeature f2old( fid2 );
397  f2old.setGeometry( g2 );
398  index.deleteFeature( f2old );
399 
400  geometries.remove( fid2 );
401 
402  if ( sanitizeDifferenceResult( g21, geometryType ) )
403  {
404  geometries.insert( fid2, g21 );
405 
406  QgsFeature f2x( fid2 );
407  f2x.setGeometry( g21 );
408  index.addFeature( f2x );
409  }
410 
411  // update our temporary copy of the geometry to what is left from it
412  g1 = g12;
413  g1engine.reset();
414  }
415 
416  ++count;
417  feedback->setProgress( count / ( double ) totalCount * 100. );
418  }
419 
420  // release some memory of structures we don't need anymore
421 
422  fids.clear();
423  index = QgsSpatialIndex();
424 
425  // load attributes
426 
427  QHash<QgsFeatureId, QgsAttributes> attributesHash;
428  it = source.getFeatures( requestOnlyAttrs );
429  while ( it.nextFeature( f ) )
430  {
431  if ( feedback->isCanceled() )
432  return;
433 
434  attributesHash.insert( f.id(), f.attributes() );
435  }
436 
437  // store stuff in the sink
438 
439  for ( auto i = geometries.constBegin(); i != geometries.constEnd(); ++i )
440  {
441  if ( feedback->isCanceled() )
442  return;
443 
444  QgsFeature outFeature( i.key() );
445  outFeature.setGeometry( i.value() );
446 
447  if ( intersectingIds.contains( i.key() ) )
448  {
449  const QList<QgsFeatureId> ids = intersectingIds.value( i.key() );
450  for ( QgsFeatureId id : ids )
451  {
452  outFeature.setAttributes( attributesHash.value( id ) );
453  sink.addFeature( outFeature, QgsFeatureSink::FastInsert );
454  }
455  }
456  else
457  {
458  outFeature.setAttributes( attributesHash.value( i.key() ) );
459  sink.addFeature( outFeature, QgsFeatureSink::FastInsert );
460  }
461  }
462 }
463 
bool addFeature(QgsFeature &feature, QgsFeatureSink::Flags flags=QgsFeatureSink::Flags()) override
Adds a feature to the index.
bool convertGeometryCollectionToSubclass(QgsWkbTypes::GeometryType geomType)
Converts geometry collection to a the desired geometry type subclass (multi-point, multi-linestring or multi-polygon).
QgsWkbTypes::Type wkbType() const SIP_HOLDGIL
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
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
QSet< QgsFeatureId > QgsFeatureIds
Definition: qgsfeatureid.h:37
QString lastError() const SIP_HOLDGIL
Returns an error string referring to the last error encountered either when this geometry was created...
Base class for providing feedback from a processing algorithm.
virtual QgsFields fields() const =0
Returns the fields associated with features in the source.
static Type multiType(Type type) SIP_HOLDGIL
Returns the multi type for a WKB type.
Definition: qgswkbtypes.h:302
QgsFeatureRequest & setInvalidGeometryCheck(InvalidGeometryCheck check)
Sets invalid geometry checking behavior.
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:62
An interface for objects which accept features via addFeature(s) methods.
qint64 QgsFeatureId
64 bit feature ids negative numbers are used for uncommitted/newly added features ...
Definition: qgsfeatureid.h:28
QgsFeatureRequest & setSubsetOfAttributes(const QgsAttributeList &attrs)
Set a subset of attributes that will be fetched.
static Type flatType(Type type) SIP_HOLDGIL
Returns the flat type for a WKB type.
Definition: qgswkbtypes.h:702
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:123
void setAttributes(const QgsAttributes &attrs)
Sets the feature&#39;s attributes.
Definition: qgsfeature.cpp:134
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:55
bool hasGeometry() const
Returns true if the feature has an associated geometry.
Definition: qgsfeature.cpp:204
int count() const
Returns number of items.
Definition: qgsfields.cpp:133
QgsFeatureRequest & setNoAttributes()
Set that no attributes will be fetched.
QgsFeatureRequest::InvalidGeometryCheck invalidGeometryCheck() const
Returns the behavior used for checking invalid geometries in input layers.
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...
const QgsAbstractGeometry * constGet() const SIP_HOLDGIL
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
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.
GeometryType
The geometry types are used to group QgsWkbTypes::Type in a coarse way.
Definition: qgswkbtypes.h:140
bool deleteFeature(const QgsFeature &feature)
Removes a feature from the index.
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:53
An interface for objects which provide features via a getFeatures method.
virtual bool addFeature(QgsFeature &feature, QgsFeatureSink::Flags flags=QgsFeatureSink::Flags())
Adds a single feature to the sink.
QgsRectangle boundingBox() const
Returns the bounding box of the geometry.
void setGeometry(const QgsGeometry &geometry)
Set the feature&#39;s geometry.
Definition: qgsfeature.cpp:144
QgsGeometry geometry
Definition: qgsfeature.h:67
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:57
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.
static GeometryType geometryType(Type type) SIP_HOLDGIL
Returns the geometry type for a WKB type, e.g., both MultiPolygon and CurvePolygon would have a Polyg...
Definition: qgswkbtypes.h:938
QgsAttributes attributes
Definition: qgsfeature.h:65
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...