QGIS API Documentation  3.21.0-Master (5b68dc587e)
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"
25 #include "qgsspatialindexutils.h"
26 
27 #include <spatialindex/SpatialIndex.h>
28 #include <QMutex>
29 #include <QMutexLocker>
30 
31 using namespace SpatialIndex;
32 
33 
34 
41 class 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 
67 class 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 
92 class 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 
162 class 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:
232  QgsFeatureIterator mFi;
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 
247 class 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 
345 QgsSpatialIndex::QgsSpatialIndex( QgsSpatialIndex::Flags flags )
346 {
347  d = new QgsSpatialIndexData( flags );
348 }
349 
350 QgsSpatialIndex::QgsSpatialIndex( const QgsFeatureIterator &fi, QgsFeedback *feedback, QgsSpatialIndex::Flags flags )
351 {
352  d = new QgsSpatialIndexData( fi, feedback, flags );
353 }
354 
356 QgsSpatialIndex::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 
362 QgsSpatialIndex::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 
383 bool 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 
393 bool 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 
407 bool 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 
426 bool 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 
492 QList<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 
505 QList<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 
531 QList<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 
562 QAtomicInt 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:205
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:124
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:834
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