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