QGIS API Documentation 3.27.0-Master (75dc696944)
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;
235 QgsSpatialIndex::Flags mFlags = QgsSpatialIndex::Flags();
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
256 QgsSpatialIndex::Flags mFlags = QgsSpatialIndex::Flags();
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#if QT_VERSION < QT_VERSION_CHECK(5, 14, 0)
333 mutable QMutex mMutex;
334#else
335 mutable QRecursiveMutex mMutex;
336#endif
337
338};
339
341
342// -------------------------------------------------------------------------
343
344
345QgsSpatialIndex::QgsSpatialIndex( QgsSpatialIndex::Flags flags )
346{
347 d = new QgsSpatialIndexData( flags );
348}
349
350QgsSpatialIndex::QgsSpatialIndex( const QgsFeatureIterator &fi, QgsFeedback *feedback, QgsSpatialIndex::Flags flags )
351{
352 d = new QgsSpatialIndexData( fi, feedback, flags );
353}
354
356QgsSpatialIndex::QgsSpatialIndex( const QgsFeatureIterator &fi, const std::function< bool( const QgsFeature & )> &callback, QgsSpatialIndex::Flags flags )
357{
358 d = new QgsSpatialIndexData( fi, nullptr, flags, &callback );
359}
361
362QgsSpatialIndex::QgsSpatialIndex( const QgsFeatureSource &source, QgsFeedback *feedback, QgsSpatialIndex::Flags flags )
363{
364 d = new QgsSpatialIndexData( source.getFeatures( QgsFeatureRequest().setNoAttributes() ), feedback, flags );
365}
366
368 : d( other.d )
369{
370}
371
373{
374}
375
377{
378 if ( this != &other )
379 d = other.d;
380 return *this;
381}
382
383bool QgsSpatialIndex::featureInfo( const QgsFeature &f, SpatialIndex::Region &r, QgsFeatureId &id )
384{
385 QgsRectangle rect;
386 if ( !featureInfo( f, rect, id ) )
387 return false;
388
390 return true;
391}
392
393bool QgsSpatialIndex::featureInfo( const QgsFeature &f, QgsRectangle &rect, QgsFeatureId &id )
394{
395 if ( !f.hasGeometry() )
396 return false;
397
398 id = f.id();
399 rect = f.geometry().boundingBox();
400
401 if ( !rect.isFinite() )
402 return false;
403
404 return true;
405}
406
407bool QgsSpatialIndex::addFeature( QgsFeature &feature, QgsFeatureSink::Flags )
408{
409 QgsRectangle rect;
410 QgsFeatureId id;
411 if ( !featureInfo( feature, rect, id ) )
412 return false;
413
414 if ( addFeature( id, rect ) )
415 {
417 {
418 const QMutexLocker locker( &d->mMutex );
419 d->mGeometries.insert( feature.id(), feature.geometry() );
420 }
421 return true;
422 }
423 return false;
424}
425
426bool QgsSpatialIndex::addFeatures( QgsFeatureList &features, QgsFeatureSink::Flags flags )
427{
428 QgsFeatureList::iterator fIt = features.begin();
429 bool result = true;
430 for ( ; fIt != features.end(); ++fIt )
431 {
432 result = result && addFeature( *fIt, flags );
433 }
434 return result;
435}
436
438{
439 QgsFeature feature( f );
440 return addFeature( feature );
441}
442
444{
445 return addFeature( id, bounds );
446}
447
449{
450 const SpatialIndex::Region r( QgsSpatialIndexUtils::rectangleToRegion( bounds ) );
451
452 const QMutexLocker locker( &d->mMutex );
453
454 // TODO: handle possible exceptions correctly
455 try
456 {
457 d->mRTree->insertData( 0, nullptr, r, FID_TO_NUMBER( id ) );
458 return true;
459 }
460 catch ( Tools::Exception &e )
461 {
462 Q_UNUSED( e )
463 QgsDebugMsg( QStringLiteral( "Tools::Exception caught: " ).arg( e.what().c_str() ) );
464 }
465 catch ( const std::exception &e )
466 {
467 Q_UNUSED( e )
468 QgsDebugMsg( QStringLiteral( "std::exception caught: " ).arg( e.what() ) );
469 }
470 catch ( ... )
471 {
472 QgsDebugMsg( QStringLiteral( "unknown spatial index exception caught" ) );
473 }
474
475 return false;
476}
477
479{
480 SpatialIndex::Region r;
481 QgsFeatureId id;
482 if ( !featureInfo( f, r, id ) )
483 return false;
484
485 const QMutexLocker locker( &d->mMutex );
486 // TODO: handle exceptions
488 d->mGeometries.remove( f.id() );
489 return d->mRTree->deleteData( r, FID_TO_NUMBER( id ) );
490}
491
492QList<QgsFeatureId> QgsSpatialIndex::intersects( const QgsRectangle &rect ) const
493{
494 QList<QgsFeatureId> list;
495 QgisVisitor visitor( list );
496
497 const SpatialIndex::Region r = QgsSpatialIndexUtils::rectangleToRegion( rect );
498
499 const QMutexLocker locker( &d->mMutex );
500 d->mRTree->intersectsWithQuery( r, visitor );
501
502 return list;
503}
504
505QList<QgsFeatureId> QgsSpatialIndex::nearestNeighbor( const QgsPointXY &point, const int neighbors, const double maxDistance ) const
506{
507 QList<QgsFeatureId> list;
508 QgisVisitor visitor( list );
509
510 double pt[2] = { point.x(), point.y() };
511 const Point p( pt, 2 );
512
513 const QMutexLocker locker( &d->mMutex );
514 QgsNearestNeighborComparator nnc( ( d->mFlags & QgsSpatialIndex::FlagStoreFeatureGeometries ) ? &d->mGeometries : nullptr,
515 point, maxDistance );
516 d->mRTree->nearestNeighborQuery( neighbors, p, visitor, nnc );
517
518 if ( maxDistance > 0 )
519 {
520 // trim features outside of max distance
521 list.erase( std::remove_if( list.begin(), list.end(),
522 [&nnc]( QgsFeatureId id )
523 {
524 return nnc.mFeaturesOutsideMaxDistance.contains( id );
525 } ), list.end() );
526 }
527
528 return list;
529}
530
531QList<QgsFeatureId> QgsSpatialIndex::nearestNeighbor( const QgsGeometry &geometry, int neighbors, double maxDistance ) const
532{
533 QList<QgsFeatureId> list;
534 QgisVisitor visitor( list );
535
536 const SpatialIndex::Region r = QgsSpatialIndexUtils::rectangleToRegion( geometry.boundingBox() );
537
538 const QMutexLocker locker( &d->mMutex );
539 QgsNearestNeighborComparator nnc( ( d->mFlags & QgsSpatialIndex::FlagStoreFeatureGeometries ) ? &d->mGeometries : nullptr,
540 geometry, maxDistance );
541 d->mRTree->nearestNeighborQuery( neighbors, r, visitor, nnc );
542
543 if ( maxDistance > 0 )
544 {
545 // trim features outside of max distance
546 list.erase( std::remove_if( list.begin(), list.end(),
547 [&nnc]( QgsFeatureId id )
548 {
549 return nnc.mFeaturesOutsideMaxDistance.contains( id );
550 } ), list.end() );
551 }
552
553 return list;
554}
555
557{
558 const QMutexLocker locker( &d->mMutex );
559 return d->mGeometries.value( id );
560}
561
562QAtomicInt QgsSpatialIndex::refs() const
563{
564 return d->ref;
565}
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).
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:45
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:125
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:59
double y
Definition: qgspointxy.h:63
Q_GADGET double x
Definition: qgspointxy.h:62
A rectangle specified with double values.
Definition: qgsrectangle.h:42
bool isFinite() const
Returns true if the rectangle has finite boundaries.
Definition: qgsrectangle.h:559
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.
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 QgsDebugMsg(str)
Definition: qgslogger.h:38