Branch data Line data Source code
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 : :
21 : : ///@cond PRIVATE
22 : :
23 : 0 : bool QgsOverlayUtils::sanitizeIntersectionResult( QgsGeometry &geom, QgsWkbTypes::GeometryType geometryType )
24 : : {
25 : 0 : 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 : 0 : 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.
34 : 0 : if ( QgsWkbTypes::flatType( geom.wkbType() ) == QgsWkbTypes::GeometryCollection )
35 : : {
36 : : // try to filter out irrelevant parts with different geometry type than what we want
37 : 0 : geom.convertGeometryCollectionToSubclass( geometryType );
38 : 0 : if ( geom.isEmpty() )
39 : 0 : return false;
40 : 0 : }
41 : :
42 : 0 : if ( QgsWkbTypes::geometryType( geom.wkbType() ) != geometryType )
43 : : {
44 : : // we can't make use of this resulting geometry
45 : 0 : 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 : 0 : geom.convertToMultiType();
51 : :
52 : 0 : return true;
53 : 0 : }
54 : :
55 : :
56 : : //! Makes sure that what came out from difference of two geometries is good to be used in the output
57 : 0 : static bool sanitizeDifferenceResult( QgsGeometry &geom, QgsWkbTypes::GeometryType geometryType )
58 : : {
59 : 0 : 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 : 0 : throw QgsProcessingException( QStringLiteral( "%1\n\n%2" ).arg( QObject::tr( "GEOS geoprocessing error: difference failed." ), geom.lastError() ) );
63 : : }
64 : :
65 : : //fix geometry collections
66 : 0 : if ( QgsWkbTypes::flatType( geom.wkbType() ) == QgsWkbTypes::GeometryCollection )
67 : : {
68 : : // try to filter out irrelevant parts with different geometry type than what we want
69 : 0 : geom.convertGeometryCollectionToSubclass( geometryType );
70 : 0 : }
71 : :
72 : :
73 : : // if geomB covers the whole source geometry, we get an empty geometry collection
74 : 0 : if ( geom.isEmpty() )
75 : 0 : 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 : 0 : geom.convertToMultiType();
80 : :
81 : 0 : return true;
82 : 0 : }
83 : :
84 : :
85 : 0 : void QgsOverlayUtils::difference( const QgsFeatureSource &sourceA, const QgsFeatureSource &sourceB, QgsFeatureSink &sink, QgsProcessingContext &context, QgsProcessingFeedback *feedback, int &count, int totalCount, QgsOverlayUtils::DifferenceOutput outputAttrs )
86 : : {
87 : 0 : QgsWkbTypes::GeometryType geometryType = QgsWkbTypes::geometryType( QgsWkbTypes::multiType( sourceA.wkbType() ) );
88 : 0 : QgsFeatureRequest requestB;
89 : 0 : requestB.setNoAttributes();
90 : 0 : if ( outputAttrs != OutputBA )
91 : 0 : requestB.setDestinationCrs( sourceA.sourceCrs(), context.transformContext() );
92 : 0 : QgsSpatialIndex indexB( sourceB.getFeatures( requestB ), feedback );
93 : :
94 : 0 : int fieldsCountA = sourceA.fields().count();
95 : 0 : int fieldsCountB = sourceB.fields().count();
96 : 0 : QgsAttributes attrs;
97 : 0 : attrs.resize( outputAttrs == OutputA ? fieldsCountA : ( fieldsCountA + fieldsCountB ) );
98 : :
99 : 0 : if ( totalCount == 0 )
100 : 0 : totalCount = 1; // avoid division by zero
101 : :
102 : 0 : QgsFeature featA;
103 : 0 : QgsFeatureRequest requestA;
104 : 0 : requestA.setInvalidGeometryCheck( context.invalidGeometryCheck() );
105 : 0 : if ( outputAttrs == OutputBA )
106 : 0 : requestA.setDestinationCrs( sourceB.sourceCrs(), context.transformContext() );
107 : 0 : QgsFeatureIterator fitA = sourceA.getFeatures( requestA );
108 : 0 : while ( fitA.nextFeature( featA ) )
109 : : {
110 : 0 : if ( feedback->isCanceled() )
111 : 0 : break;
112 : :
113 : 0 : if ( featA.hasGeometry() )
114 : : {
115 : 0 : QgsGeometry geom( featA.geometry() );
116 : 0 : QgsFeatureIds intersects = qgis::listToSet( indexB.intersects( geom.boundingBox() ) );
117 : :
118 : 0 : QgsFeatureRequest request;
119 : 0 : request.setFilterFids( intersects );
120 : 0 : request.setNoAttributes();
121 : 0 : if ( outputAttrs != OutputBA )
122 : 0 : request.setDestinationCrs( sourceA.sourceCrs(), context.transformContext() );
123 : :
124 : 0 : std::unique_ptr< QgsGeometryEngine > engine;
125 : 0 : if ( !intersects.isEmpty() )
126 : : {
127 : : // use prepared geometries for faster intersection tests
128 : 0 : engine.reset( QgsGeometry::createGeometryEngine( geom.constGet() ) );
129 : 0 : engine->prepareGeometry();
130 : 0 : }
131 : :
132 : 0 : QVector<QgsGeometry> geometriesB;
133 : 0 : QgsFeature featB;
134 : 0 : QgsFeatureIterator fitB = sourceB.getFeatures( request );
135 : 0 : while ( fitB.nextFeature( featB ) )
136 : : {
137 : 0 : if ( feedback->isCanceled() )
138 : 0 : break;
139 : :
140 : 0 : if ( engine->intersects( featB.geometry().constGet() ) )
141 : 0 : geometriesB << featB.geometry();
142 : : }
143 : :
144 : 0 : if ( !geometriesB.isEmpty() )
145 : : {
146 : 0 : QgsGeometry geomB = QgsGeometry::unaryUnion( geometriesB );
147 : 0 : 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 : 0 : throw QgsProcessingException( QStringLiteral( "%1\n\n%2" ).arg( QObject::tr( "GEOS geoprocessing error: unary union failed." ), geomB.lastError() ) );
155 : : }
156 : 0 : geom = geom.difference( geomB );
157 : 0 : }
158 : :
159 : 0 : if ( !sanitizeDifferenceResult( geom, geometryType ) )
160 : 0 : continue;
161 : :
162 : 0 : const QgsAttributes attrsA( featA.attributes() );
163 : 0 : switch ( outputAttrs )
164 : : {
165 : : case OutputA:
166 : 0 : attrs = attrsA;
167 : 0 : break;
168 : : case OutputAB:
169 : 0 : for ( int i = 0; i < fieldsCountA; ++i )
170 : 0 : attrs[i] = attrsA[i];
171 : 0 : break;
172 : : case OutputBA:
173 : 0 : for ( int i = 0; i < fieldsCountA; ++i )
174 : 0 : attrs[i + fieldsCountB] = attrsA[i];
175 : 0 : break;
176 : : }
177 : :
178 : 0 : QgsFeature outFeat;
179 : 0 : outFeat.setGeometry( geom );
180 : 0 : outFeat.setAttributes( attrs );
181 : 0 : sink.addFeature( outFeat, QgsFeatureSink::FastInsert );
182 : 0 : }
183 : : else
184 : : {
185 : : // TODO: should we write out features that do not have geometry?
186 : 0 : sink.addFeature( featA, QgsFeatureSink::FastInsert );
187 : : }
188 : :
189 : 0 : ++count;
190 : 0 : feedback->setProgress( count / ( double ) totalCount * 100. );
191 : : }
192 : 0 : }
193 : :
194 : :
195 : 0 : 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 : : {
197 : 0 : QgsWkbTypes::GeometryType geometryType = QgsWkbTypes::geometryType( QgsWkbTypes::multiType( sourceA.wkbType() ) );
198 : 0 : int attrCount = fieldIndicesA.count() + fieldIndicesB.count();
199 : :
200 : 0 : QgsFeatureRequest request;
201 : 0 : request.setNoAttributes();
202 : 0 : request.setDestinationCrs( sourceA.sourceCrs(), context.transformContext() );
203 : :
204 : 0 : QgsFeature outFeat;
205 : 0 : QgsSpatialIndex indexB( sourceB.getFeatures( request ), feedback );
206 : :
207 : 0 : if ( totalCount == 0 )
208 : 0 : totalCount = 1; // avoid division by zero
209 : :
210 : 0 : QgsFeature featA;
211 : 0 : QgsFeatureIterator fitA = sourceA.getFeatures( QgsFeatureRequest().setSubsetOfAttributes( fieldIndicesA ) );
212 : 0 : while ( fitA.nextFeature( featA ) )
213 : : {
214 : 0 : if ( feedback->isCanceled() )
215 : 0 : break;
216 : :
217 : 0 : if ( !featA.hasGeometry() )
218 : 0 : continue;
219 : :
220 : 0 : QgsGeometry geom( featA.geometry() );
221 : 0 : QgsFeatureIds intersects = qgis::listToSet( indexB.intersects( geom.boundingBox() ) );
222 : :
223 : 0 : QgsFeatureRequest request;
224 : 0 : request.setFilterFids( intersects );
225 : 0 : request.setDestinationCrs( sourceA.sourceCrs(), context.transformContext() );
226 : 0 : request.setSubsetOfAttributes( fieldIndicesB );
227 : :
228 : 0 : std::unique_ptr< QgsGeometryEngine > engine;
229 : 0 : if ( !intersects.isEmpty() )
230 : : {
231 : : // use prepared geometries for faster intersection tests
232 : 0 : engine.reset( QgsGeometry::createGeometryEngine( geom.constGet() ) );
233 : 0 : engine->prepareGeometry();
234 : 0 : }
235 : :
236 : 0 : QgsAttributes outAttributes( attrCount );
237 : 0 : const QgsAttributes attrsA( featA.attributes() );
238 : 0 : for ( int i = 0; i < fieldIndicesA.count(); ++i )
239 : 0 : outAttributes[i] = attrsA[fieldIndicesA[i]];
240 : :
241 : 0 : QgsFeature featB;
242 : 0 : QgsFeatureIterator fitB = sourceB.getFeatures( request );
243 : 0 : while ( fitB.nextFeature( featB ) )
244 : : {
245 : 0 : if ( feedback->isCanceled() )
246 : 0 : break;
247 : :
248 : 0 : QgsGeometry tmpGeom( featB.geometry() );
249 : 0 : if ( !engine->intersects( tmpGeom.constGet() ) )
250 : 0 : continue;
251 : :
252 : 0 : QgsGeometry intGeom = geom.intersection( tmpGeom );
253 : 0 : if ( !sanitizeIntersectionResult( intGeom, geometryType ) )
254 : 0 : continue;
255 : :
256 : 0 : const QgsAttributes attrsB( featB.attributes() );
257 : 0 : for ( int i = 0; i < fieldIndicesB.count(); ++i )
258 : 0 : outAttributes[fieldIndicesA.count() + i] = attrsB[fieldIndicesB[i]];
259 : :
260 : 0 : outFeat.setGeometry( intGeom );
261 : 0 : outFeat.setAttributes( outAttributes );
262 : 0 : sink.addFeature( outFeat, QgsFeatureSink::FastInsert );
263 : 0 : }
264 : :
265 : 0 : ++count;
266 : 0 : feedback->setProgress( count / ( double ) totalCount * 100. );
267 : 0 : }
268 : 0 : }
269 : :
270 : 0 : void QgsOverlayUtils::resolveOverlaps( const QgsFeatureSource &source, QgsFeatureSink &sink, QgsProcessingFeedback *feedback )
271 : : {
272 : 0 : int count = 0;
273 : 0 : int totalCount = source.featureCount();
274 : 0 : if ( totalCount == 0 )
275 : 0 : return; // nothing to do here
276 : :
277 : 0 : QgsFeatureId newFid = -1;
278 : :
279 : 0 : QgsWkbTypes::GeometryType geometryType = QgsWkbTypes::geometryType( QgsWkbTypes::multiType( source.wkbType() ) );
280 : :
281 : 0 : QgsFeatureRequest requestOnlyGeoms;
282 : 0 : requestOnlyGeoms.setNoAttributes();
283 : :
284 : 0 : QgsFeatureRequest requestOnlyAttrs;
285 : 0 : requestOnlyAttrs.setFlags( QgsFeatureRequest::NoGeometry );
286 : :
287 : 0 : QgsFeatureRequest requestOnlyIds;
288 : 0 : requestOnlyIds.setFlags( QgsFeatureRequest::NoGeometry );
289 : 0 : 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 : 0 : QgsFeature f;
293 : 0 : QSet<QgsFeatureId> fids;
294 : 0 : QgsFeatureIterator it = source.getFeatures( requestOnlyIds );
295 : 0 : while ( it.nextFeature( f ) )
296 : : {
297 : 0 : if ( feedback->isCanceled() )
298 : 0 : return;
299 : :
300 : 0 : fids.insert( f.id() );
301 : : }
302 : :
303 : 0 : QHash<QgsFeatureId, QgsGeometry> geometries;
304 : 0 : QgsSpatialIndex index;
305 : 0 : QHash<QgsFeatureId, QList<QgsFeatureId> > intersectingIds; // which features overlap a particular area
306 : :
307 : : // resolve intersections
308 : :
309 : 0 : it = source.getFeatures( requestOnlyGeoms );
310 : 0 : while ( it.nextFeature( f ) )
311 : : {
312 : 0 : if ( feedback->isCanceled() )
313 : 0 : return;
314 : :
315 : 0 : QgsFeatureId fid1 = f.id();
316 : 0 : QgsGeometry g1 = f.geometry();
317 : 0 : std::unique_ptr< QgsGeometryEngine > g1engine;
318 : :
319 : 0 : geometries.insert( fid1, g1 );
320 : 0 : index.addFeature( f );
321 : :
322 : 0 : QgsRectangle bbox( f.geometry().boundingBox() );
323 : 0 : const QList<QgsFeatureId> ids = index.intersects( bbox );
324 : 0 : for ( QgsFeatureId fid2 : ids )
325 : : {
326 : 0 : if ( fid1 == fid2 )
327 : 0 : continue;
328 : :
329 : 0 : if ( !g1engine )
330 : : {
331 : : // use prepared geometries for faster intersection tests
332 : 0 : g1engine.reset( QgsGeometry::createGeometryEngine( g1.constGet() ) );
333 : 0 : g1engine->prepareGeometry();
334 : 0 : }
335 : :
336 : 0 : QgsGeometry g2 = geometries.value( fid2 );
337 : 0 : if ( !g1engine->intersects( g2.constGet() ) )
338 : 0 : continue;
339 : :
340 : 0 : QgsGeometry geomIntersection = g1.intersection( g2 );
341 : 0 : if ( !sanitizeIntersectionResult( geomIntersection, geometryType ) )
342 : 0 : continue;
343 : :
344 : : //
345 : : // add intersection geometry
346 : : //
347 : :
348 : : // figure out new fid
349 : 0 : while ( fids.contains( newFid ) )
350 : 0 : --newFid;
351 : 0 : fids.insert( newFid );
352 : :
353 : 0 : geometries.insert( newFid, geomIntersection );
354 : 0 : QgsFeature fx( newFid );
355 : 0 : fx.setGeometry( geomIntersection );
356 : :
357 : 0 : 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 : 0 : QList<QgsFeatureId> lst;
362 : 0 : if ( intersectingIds.contains( fid1 ) )
363 : 0 : lst << intersectingIds.value( fid1 );
364 : : else
365 : 0 : lst << fid1;
366 : 0 : if ( intersectingIds.contains( fid2 ) )
367 : 0 : lst << intersectingIds.value( fid2 );
368 : : else
369 : 0 : lst << fid2;
370 : 0 : intersectingIds.insert( newFid, lst );
371 : :
372 : : //
373 : : // update f1
374 : : //
375 : :
376 : 0 : QgsGeometry g12 = g1.difference( g2 );
377 : :
378 : 0 : index.deleteFeature( f );
379 : 0 : geometries.remove( fid1 );
380 : :
381 : 0 : if ( sanitizeDifferenceResult( g12, geometryType ) )
382 : : {
383 : 0 : geometries.insert( fid1, g12 );
384 : :
385 : 0 : QgsFeature f1x( fid1 );
386 : 0 : f1x.setGeometry( g12 );
387 : 0 : index.addFeature( f1x );
388 : 0 : }
389 : :
390 : : //
391 : : // update f2
392 : : //
393 : :
394 : 0 : QgsGeometry g21 = g2.difference( g1 );
395 : :
396 : 0 : QgsFeature f2old( fid2 );
397 : 0 : f2old.setGeometry( g2 );
398 : 0 : index.deleteFeature( f2old );
399 : :
400 : 0 : geometries.remove( fid2 );
401 : :
402 : 0 : if ( sanitizeDifferenceResult( g21, geometryType ) )
403 : : {
404 : 0 : geometries.insert( fid2, g21 );
405 : :
406 : 0 : QgsFeature f2x( fid2 );
407 : 0 : f2x.setGeometry( g21 );
408 : 0 : index.addFeature( f2x );
409 : 0 : }
410 : :
411 : : // update our temporary copy of the geometry to what is left from it
412 : 0 : g1 = g12;
413 : 0 : g1engine.reset();
414 : 0 : }
415 : :
416 : 0 : ++count;
417 : 0 : feedback->setProgress( count / ( double ) totalCount * 100. );
418 : 0 : }
419 : :
420 : : // release some memory of structures we don't need anymore
421 : :
422 : 0 : fids.clear();
423 : 0 : index = QgsSpatialIndex();
424 : :
425 : : // load attributes
426 : :
427 : 0 : QHash<QgsFeatureId, QgsAttributes> attributesHash;
428 : 0 : it = source.getFeatures( requestOnlyAttrs );
429 : 0 : while ( it.nextFeature( f ) )
430 : : {
431 : 0 : if ( feedback->isCanceled() )
432 : 0 : return;
433 : :
434 : 0 : attributesHash.insert( f.id(), f.attributes() );
435 : : }
436 : :
437 : : // store stuff in the sink
438 : :
439 : 0 : for ( auto i = geometries.constBegin(); i != geometries.constEnd(); ++i )
440 : : {
441 : 0 : if ( feedback->isCanceled() )
442 : 0 : return;
443 : :
444 : 0 : QgsFeature outFeature( i.key() );
445 : 0 : outFeature.setGeometry( i.value() );
446 : :
447 : 0 : if ( intersectingIds.contains( i.key() ) )
448 : : {
449 : 0 : const QList<QgsFeatureId> ids = intersectingIds.value( i.key() );
450 : 0 : for ( QgsFeatureId id : ids )
451 : : {
452 : 0 : outFeature.setAttributes( attributesHash.value( id ) );
453 : 0 : sink.addFeature( outFeature, QgsFeatureSink::FastInsert );
454 : : }
455 : 0 : }
456 : : else
457 : : {
458 : 0 : outFeature.setAttributes( attributesHash.value( i.key() ) );
459 : 0 : sink.addFeature( outFeature, QgsFeatureSink::FastInsert );
460 : : }
461 : 0 : }
462 : 0 : }
463 : :
464 : : ///@endcond PRIVATE
|
