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