2

To be more specific: I have a feature class for which some features have an entry for a variable ('KnownBj') and some not (have a 9999 entry). I want for each feature to find the 5 nearest neighbors that have an entry in 'KnownBj'. What I currently do is to prepare two dictionaries - A has all the features and B only the features with a 'KnownBj' .I add B to a SpatialIndex, and loop through A searching for neighbours with .NearestNeighbor:

canvas = qgis.utils.iface.mapCanvas()
layer = canvas.currentLayer()
iter = layer.selectedFeatures()
feature_dict = {f.id(): f for f in iter}
featDictKnownBj = {k: v for k, v in feature_dict.iteritems() if v[KnownBj] != '9999' }

#Prepare Spatial Index
SpatIndex = QgsSpatialIndex()
map(SpatIndex.insertFeature, featDictKnownBj.values())

#SEARCH NEIGHBOURS
NNDict = {}
for f in feature_dict.values():
    NNs = SpatIndex.nearestNeighbor(f.geometry().centroid().asPoint(),5)

Now this seems to work, my question is - am I missing sth? Should I add the current feature to the Spatial Index, then search, then remove it and so for all, or should this be perfectly ok?

IvDoch
  • 33
  • 2
  • Spatial index just can make things faster. It stores the bounding boxes of features. Bboxes are fast to query but they are not accurate and the final query will use original geometries. – user30184 Oct 09 '16 at 11:18

1 Answers1

1

As user30184 comments, in order to find the five true nearest neighbours you need to consider the original geometries, and the five true nearest neighbours do not need to be among the five features returned by QgsSpatialIndex.nearestNeighbor.

To get the five true nearest neighbours you could go through this process:

  1. Find the distances to the geometries of the five features returned by QgsSpatialIndex.nearestNeighbor using using QgsGeometry.distance.
  2. Use the largest of these distances to create a rectangle that is used to select all possible candidates using QgsSpatialIndex.intersects.
  3. Find the five closest of these candidate features using QgsGeometry.distance.

In https://gis.stackexchange.com/a/257408 there is code for a similar case.

Håvard Tveite
  • 3,256
  • 17
  • 32