QGIS API Documentation 3.37.0-Master (fdefdf9c27f)
qgsspatialindex.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsspatialindex.cpp - wrapper class for spatial index library
3 ----------------------
4 begin : December 2006
5 copyright : (C) 2006 by Martin Dobias
6 email : wonder.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 "qgsspatialindex.h"
17
18#include "qgsgeometry.h"
19#include "qgsfeature.h"
20#include "qgsfeatureiterator.h"
21#include "qgsrectangle.h"
22#include "qgslogger.h"
23#include "qgsfeaturesource.h"
24#include "qgsfeedback.h"
26
27#include <spatialindex/SpatialIndex.h>
28#include <QMutex>
29#include <QMutexLocker>
30
31using namespace SpatialIndex;
32
33
34
41class QgisVisitor : public SpatialIndex::IVisitor
42{
43 public:
44 explicit QgisVisitor( QList<QgsFeatureId> &list )
45 : mList( list ) {}
46
47 void visitNode( const INode &n ) override
48 { Q_UNUSED( n ) }
49
50 void visitData( const IData &d ) override
51 {
52 mList.append( d.getIdentifier() );
53 }
54
55 void visitData( std::vector<const IData *> &v ) override
56 { Q_UNUSED( v ) }
57
58 private:
59 QList<QgsFeatureId> &mList;
60};
61
67class QgsSpatialIndexCopyVisitor : public SpatialIndex::IVisitor
68{
69 public:
70 explicit QgsSpatialIndexCopyVisitor( SpatialIndex::ISpatialIndex *newIndex )
71 : mNewIndex( newIndex ) {}
72
73 void visitNode( const INode &n ) override
74 { Q_UNUSED( n ) }
75
76 void visitData( const IData &d ) override
77 {
78 SpatialIndex::IShape *shape = nullptr;
79 d.getShape( &shape );
80 mNewIndex->insertData( 0, nullptr, *shape, d.getIdentifier() );
81 delete shape;
82 }
83
84 void visitData( std::vector<const IData *> &v ) override
85 { Q_UNUSED( v ) }
86
87 private:
88 SpatialIndex::ISpatialIndex *mNewIndex = nullptr;
89};
90
92class QgsNearestNeighborComparator : public INearestNeighborComparator
93{
94 public:
95
96 QgsNearestNeighborComparator( const QHash< QgsFeatureId, QgsGeometry > *geometries, const QgsPointXY &point, double maxDistance )
97 : mGeometries( geometries )
98 , mGeom( QgsGeometry::fromPointXY( point ) )
99 , mMaxDistance( maxDistance )
100 {
101 }
102
103 QgsNearestNeighborComparator( const QHash< QgsFeatureId, QgsGeometry > *geometries, const QgsGeometry &geometry, double maxDistance )
104 : mGeometries( geometries )
105 , mGeom( geometry )
106 , mMaxDistance( maxDistance )
107 {
108 }
109
110 const QHash< QgsFeatureId, QgsGeometry > *mGeometries = nullptr;
111 QgsGeometry mGeom;
112 double mMaxDistance = 0;
113 QSet< QgsFeatureId > mFeaturesOutsideMaxDistance;
114
115 double getMinimumDistance( const IShape &query, const IShape &entry ) override
116 {
117 return query.getMinimumDistance( entry );
118 }
119
120 double getMinimumDistance( const IShape &query, const IData &data ) override
121 {
122 // start with the default implementation, which gives distance to bounding box only
123 IShape *pS;
124 data.getShape( &pS );
125 double dist = query.getMinimumDistance( *pS );
126 delete pS;
127
128 // if doing exact distance search, AND either no max distance specified OR the
129 // distance to the bounding box is less than the max distance, calculate the exact
130 // distance now.
131 // (note: if bounding box is already greater than the distance away then max distance, there's no
132 // point doing this more expensive calculation, since we can't possibly use this feature anyway!)
133 if ( mGeometries && ( mMaxDistance <= 0.0 || dist <= mMaxDistance ) )
134 {
135 const QgsGeometry other = mGeometries->value( data.getIdentifier() );
136 dist = other.distance( mGeom );
137 }
138
139 if ( mMaxDistance > 0 && dist > mMaxDistance )
140 {
141 // feature is outside of maximum distance threshold. Flag it,
142 // but "trick" libspatialindex into considering it as just outside
143 // our max distance region. This means if there's no other closer features (i.e.,
144 // within our actual maximum search distance), the libspatialindex
145 // nearest neighbor test will use this feature and not needlessly continue hunting
146 // through the remaining more distant features in the index.
147 // TODO: add proper API to libspatialindex to allow a maximum search distance in
148 // nearest neighbor tests
149 mFeaturesOutsideMaxDistance.insert( data.getIdentifier() );
150 return mMaxDistance + 0.00000001;
151 }
152 return dist;
153 }
154};
155
162class QgsFeatureIteratorDataStream : public IDataStream
163{
164 public:
166 explicit QgsFeatureIteratorDataStream( const QgsFeatureIterator &fi, QgsFeedback *feedback = nullptr, QgsSpatialIndex::Flags flags = QgsSpatialIndex::Flags(),
167 const std::function< bool( const QgsFeature & ) > *callback = nullptr )
168 : mFi( fi )
169 , mFeedback( feedback )
170 , mFlags( flags )
171 , mCallback( callback )
172 {
173 readNextEntry();
174 }
175
176 ~QgsFeatureIteratorDataStream() override
177 {
178 delete mNextData;
179 }
180
182 IData *getNext() override
183 {
184 if ( mFeedback && mFeedback->isCanceled() )
185 return nullptr;
186
187 RTree::Data *ret = mNextData;
188 mNextData = nullptr;
189 readNextEntry();
190 return ret;
191 }
192
194 bool hasNext() override { return nullptr != mNextData; }
195
197 uint32_t size() override { Q_ASSERT( false && "not available" ); return 0; }
198
200 void rewind() override { Q_ASSERT( false && "not available" ); }
201
202 QHash< QgsFeatureId, QgsGeometry > geometries;
203
204 protected:
205 void readNextEntry()
206 {
207 QgsFeature f;
208 SpatialIndex::Region r;
209 QgsFeatureId id;
210 while ( mFi.nextFeature( f ) )
211 {
212 if ( mCallback )
213 {
214 const bool res = ( *mCallback )( f );
215 if ( !res )
216 {
217 mNextData = nullptr;
218 return;
219 }
220 }
221 if ( QgsSpatialIndex::featureInfo( f, r, id ) )
222 {
223 mNextData = new RTree::Data( 0, nullptr, r, id );
225 geometries.insert( f.id(), f.geometry() );
226 return;
227 }
228 }
229 }
230
231 private:
233 RTree::Data *mNextData = nullptr;
234 QgsFeedback *mFeedback = nullptr;
236 const std::function< bool( const QgsFeature & ) > *mCallback = nullptr;
237
238};
239
240
247class QgsSpatialIndexData : public QSharedData
248{
249 public:
250 QgsSpatialIndexData( QgsSpatialIndex::Flags flags )
251 : mFlags( flags )
252 {
253 initTree();
254 }
255
257
258 QHash< QgsFeatureId, QgsGeometry > mGeometries;
259
268 explicit QgsSpatialIndexData( const QgsFeatureIterator &fi, QgsFeedback *feedback = nullptr, QgsSpatialIndex::Flags flags = QgsSpatialIndex::Flags(),
269 const std::function< bool( const QgsFeature & ) > *callback = nullptr )
270 : mFlags( flags )
271 {
272 QgsFeatureIteratorDataStream fids( fi, feedback, mFlags, callback );
273 initTree( &fids );
275 mGeometries = fids.geometries;
276 }
277
278 QgsSpatialIndexData( const QgsSpatialIndexData &other )
279 : QSharedData( other )
280 , mFlags( other.mFlags )
281 , mGeometries( other.mGeometries )
282 {
283 const QMutexLocker locker( &other.mMutex );
284
285 initTree();
286
287 // copy R-tree data one by one (is there a faster way??)
288 double low[] = { std::numeric_limits<double>::lowest(), std::numeric_limits<double>::lowest() };
289 double high[] = { std::numeric_limits<double>::max(), std::numeric_limits<double>::max() };
290 const SpatialIndex::Region query( low, high, 2 );
291 QgsSpatialIndexCopyVisitor visitor( mRTree );
292 other.mRTree->intersectsWithQuery( query, visitor );
293 }
294
295 ~QgsSpatialIndexData()
296 {
297 delete mRTree;
298 delete mStorage;
299 }
300
301 QgsSpatialIndexData &operator=( const QgsSpatialIndexData &rh ) = delete;
302
303 void initTree( IDataStream *inputStream = nullptr )
304 {
305 // for now only memory manager
306 mStorage = StorageManager::createNewMemoryStorageManager();
307
308 // R-Tree parameters
309 const double fillFactor = 0.7;
310 const unsigned long indexCapacity = 10;
311 const unsigned long leafCapacity = 10;
312 const unsigned long dimension = 2;
313 const RTree::RTreeVariant variant = RTree::RV_RSTAR;
314
315 // create R-tree
316 SpatialIndex::id_type indexId;
317
318 if ( inputStream && inputStream->hasNext() )
319 mRTree = RTree::createAndBulkLoadNewRTree( RTree::BLM_STR, *inputStream, *mStorage, fillFactor, indexCapacity,
320 leafCapacity, dimension, variant, indexId );
321 else
322 mRTree = RTree::createNewRTree( *mStorage, fillFactor, indexCapacity,
323 leafCapacity, dimension, variant, indexId );
324 }
325
327 SpatialIndex::IStorageManager *mStorage = nullptr;
328
330 SpatialIndex::ISpatialIndex *mRTree = nullptr;
331
332 mutable QRecursiveMutex mMutex;
333
334};
335
337
338// -------------------------------------------------------------------------
339
340
342{
343 d = new QgsSpatialIndexData( flags );
344}
345
347{
348 d = new QgsSpatialIndexData( fi, feedback, flags );
349}
350
352QgsSpatialIndex::QgsSpatialIndex( const QgsFeatureIterator &fi, const std::function< bool( const QgsFeature & )> &callback, QgsSpatialIndex::Flags flags )
353{
354 d = new QgsSpatialIndexData( fi, nullptr, flags, &callback );
355}
357
359{
360 d = new QgsSpatialIndexData( source.getFeatures( QgsFeatureRequest().setNoAttributes() ), feedback, flags );
361}
362
364 : d( other.d )
365{
366}
367
369{
370}
371
373{
374 if ( this != &other )
375 d = other.d;
376 return *this;
377}
378
379bool QgsSpatialIndex::featureInfo( const QgsFeature &f, SpatialIndex::Region &r, QgsFeatureId &id )
380{
381 QgsRectangle rect;
382 if ( !featureInfo( f, rect, id ) )
383 return false;
384
386 return true;
387}
388
389bool QgsSpatialIndex::featureInfo( const QgsFeature &f, QgsRectangle &rect, QgsFeatureId &id )
390{
391 if ( !f.hasGeometry() )
392 return false;
393
394 id = f.id();
395 rect = f.geometry().boundingBox();
396
397 if ( !rect.isFinite() )
398 return false;
399
400 return true;
401}
402
404{
405 QgsRectangle rect;
406 QgsFeatureId id;
407 if ( !featureInfo( feature, rect, id ) )
408 return false;
409
410 if ( addFeature( id, rect ) )
411 {
413 {
414 const QMutexLocker locker( &d->mMutex );
415 d->mGeometries.insert( feature.id(), feature.geometry() );
416 }
417 return true;
418 }
419 return false;
420}
421
423{
424 QgsFeatureList::iterator fIt = features.begin();
425 bool result = true;
426 for ( ; fIt != features.end(); ++fIt )
427 {
428 result = result && addFeature( *fIt, flags );
429 }
430 return result;
431}
432
434{
435 QgsFeature feature( f );
436 return addFeature( feature );
437}
438
440{
441 return addFeature( id, bounds );
442}
443
445{
446 const SpatialIndex::Region r( QgsSpatialIndexUtils::rectangleToRegion( bounds ) );
447
448 const QMutexLocker locker( &d->mMutex );
449
450 // TODO: handle possible exceptions correctly
451 try
452 {
453 d->mRTree->insertData( 0, nullptr, r, FID_TO_NUMBER( id ) );
454 return true;
455 }
456 catch ( Tools::Exception &e )
457 {
458 Q_UNUSED( e )
459 QgsDebugError( QStringLiteral( "Tools::Exception caught: " ).arg( e.what().c_str() ) );
460 }
461 catch ( const std::exception &e )
462 {
463 Q_UNUSED( e )
464 QgsDebugError( QStringLiteral( "std::exception caught: " ).arg( e.what() ) );
465 }
466 catch ( ... )
467 {
468 QgsDebugError( QStringLiteral( "unknown spatial index exception caught" ) );
469 }
470
471 return false;
472}
473
475{
476 SpatialIndex::Region r;
477 QgsFeatureId id;
478 if ( !featureInfo( f, r, id ) )
479 return false;
480
481 const QMutexLocker locker( &d->mMutex );
482 // TODO: handle exceptions
484 d->mGeometries.remove( f.id() );
485 return d->mRTree->deleteData( r, FID_TO_NUMBER( id ) );
486}
487
489{
490 const SpatialIndex::Region r = QgsSpatialIndexUtils::rectangleToRegion( bounds );
491
492 const QMutexLocker locker( &d->mMutex );
493 // TODO: handle exceptions (if the comment in the other ::deleteFeature implementation is correct!)
495 d->mGeometries.remove( id );
496 return d->mRTree->deleteData( r, FID_TO_NUMBER( id ) );
497}
498
499QList<QgsFeatureId> QgsSpatialIndex::intersects( const QgsRectangle &rect ) const
500{
501 QList<QgsFeatureId> list;
502 QgisVisitor visitor( list );
503
504 const SpatialIndex::Region r = QgsSpatialIndexUtils::rectangleToRegion( rect );
505
506 const QMutexLocker locker( &d->mMutex );
507 d->mRTree->intersectsWithQuery( r, visitor );
508
509 return list;
510}
511
512QList<QgsFeatureId> QgsSpatialIndex::nearestNeighbor( const QgsPointXY &point, const int neighbors, const double maxDistance ) const
513{
514 QList<QgsFeatureId> list;
515 QgisVisitor visitor( list );
516
517 double pt[2] = { point.x(), point.y() };
518 const Point p( pt, 2 );
519
520 const QMutexLocker locker( &d->mMutex );
521 QgsNearestNeighborComparator nnc( ( d->mFlags & QgsSpatialIndex::FlagStoreFeatureGeometries ) ? &d->mGeometries : nullptr,
522 point, maxDistance );
523 d->mRTree->nearestNeighborQuery( neighbors, p, visitor, nnc );
524
525 if ( maxDistance > 0 )
526 {
527 // trim features outside of max distance
528 list.erase( std::remove_if( list.begin(), list.end(),
529 [&nnc]( QgsFeatureId id )
530 {
531 return nnc.mFeaturesOutsideMaxDistance.contains( id );
532 } ), list.end() );
533 }
534
535 return list;
536}
537
538QList<QgsFeatureId> QgsSpatialIndex::nearestNeighbor( const QgsGeometry &geometry, int neighbors, double maxDistance ) const
539{
540 QList<QgsFeatureId> list;
541 QgisVisitor visitor( list );
542
543 const SpatialIndex::Region r = QgsSpatialIndexUtils::rectangleToRegion( geometry.boundingBox() );
544
545 const QMutexLocker locker( &d->mMutex );
546 QgsNearestNeighborComparator nnc( ( d->mFlags & QgsSpatialIndex::FlagStoreFeatureGeometries ) ? &d->mGeometries : nullptr,
547 geometry, maxDistance );
548 d->mRTree->nearestNeighborQuery( neighbors, r, visitor, nnc );
549
550 if ( maxDistance > 0 )
551 {
552 // trim features outside of max distance
553 list.erase( std::remove_if( list.begin(), list.end(),
554 [&nnc]( QgsFeatureId id )
555 {
556 return nnc.mFeaturesOutsideMaxDistance.contains( id );
557 } ), list.end() );
558 }
559
560 return list;
561}
562
564{
565 const QMutexLocker locker( &d->mMutex );
566 return d->mGeometries.value( id );
567}
568
569QAtomicInt QgsSpatialIndex::refs() const
570{
571 return d->ref;
572}
Custom visitor that adds found features to list.
void visitNode(const INode &n) override
void visitData(std::vector< const IData * > &v) override
QgisVisitor(QList< QgsFeatureId > &list)
void visitData(const IData &d) override
Wrapper for iterator of features from vector data provider or vector layer.
This class wraps a request for features to a vector layer (or directly its vector data provider).
QFlags< Flag > Flags
An interface for objects which provide features via a getFeatures method.
virtual QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest()) const =0
Returns an iterator for the features in the source.
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition: qgsfeature.h:56
QgsGeometry geometry
Definition: qgsfeature.h:67
bool hasGeometry() const
Returns true if the feature has an associated geometry.
Definition: qgsfeature.cpp:230
Q_GADGET QgsFeatureId id
Definition: qgsfeature.h:64
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition: qgsfeedback.h:44
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:162
double distance(const QgsGeometry &geom) const
Returns the minimum distance between this geometry and another geometry.
QgsRectangle boundingBox() const
Returns the bounding box of the geometry.
A class to represent a 2D point.
Definition: qgspointxy.h:60
double y
Definition: qgspointxy.h:64
Q_GADGET double x
Definition: qgspointxy.h:63
A rectangle specified with double values.
Definition: qgsrectangle.h:42
bool isFinite() const
Returns true if the rectangle has finite boundaries.
Definition: qgsrectangle.h:588
void visitData(std::vector< const IData * > &v) override
void visitData(const IData &d) override
QgsSpatialIndexCopyVisitor(SpatialIndex::ISpatialIndex *newIndex)
void visitNode(const INode &n) override
static SpatialIndex::Region rectangleToRegion(const QgsRectangle &rectangle)
Converts a QGIS rectangle to a SpatialIndex region.
A spatial index for QgsFeature objects.
bool addFeatures(QgsFeatureList &features, QgsFeatureSink::Flags flags=QgsFeatureSink::Flags()) override
Adds a list of features to the index.
@ FlagStoreFeatureGeometries
Indicates that the spatial index should also store feature geometries. This requires more memory,...
QgsSpatialIndex & operator=(const QgsSpatialIndex &other)
Implement assignment operator.
QgsSpatialIndex(QgsSpatialIndex::Flags flags=QgsSpatialIndex::Flags())
Constructor for QgsSpatialIndex.
QAtomicInt refs() const
Gets reference count - just for debugging!
QList< QgsFeatureId > nearestNeighbor(const QgsPointXY &point, int neighbors=1, double maxDistance=0) const
Returns nearest neighbors to a point.
QList< QgsFeatureId > intersects(const QgsRectangle &rectangle) const
Returns a list of features with a bounding box which intersects the specified rectangle.
~QgsSpatialIndex() override
Destructor finalizes work with spatial index.
bool addFeature(QgsFeature &feature, QgsFeatureSink::Flags flags=QgsFeatureSink::Flags()) override
Adds a feature to the index.
QFlags< Flag > Flags
Q_DECL_DEPRECATED bool insertFeature(const QgsFeature &feature)
Adds a feature to the index.
QgsGeometry geometry(QgsFeatureId id) const
Returns the stored geometry for the indexed feature with matching id.
bool deleteFeature(const QgsFeature &feature)
Removes a feature from the index.
QList< QgsFeature > QgsFeatureList
Definition: qgsfeature.h:917
qint64 QgsFeatureId
64 bit feature ids negative numbers are used for uncommitted/newly added features
Definition: qgsfeatureid.h:28
#define FID_TO_NUMBER(fid)
Definition: qgsfeatureid.h:32
#define QgsDebugError(str)
Definition: qgslogger.h:38