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