Does QgsSpatialIndex() support Z values, respectively 3D distances?

Question

I am wondering if QgsSpatialIndex() does support Z values, respectively 3D distances?

If too short to read:

Looking at the underlying libspatialindex I cannot find information neither. However, rtree does support 3D:

As of Rtree version 0.5.0, you can create 3D (actually kD) indexes. The following is a 3D index that is to be stored on disk. Persisted indexes are stored on disk using two files – an index file (.idx) and a data (.dat) file. You can modify the extensions these files use by altering the properties of the index at instantiation time. The following creates a 3D index that is stored on disk as the files 3d_index.data and 3d_index.index

Back to QgsSpatialIndex(), which

Creates an empty R-tree index.

let's take .nearestNeighbor() for example to find nearest 3D points:

The docs state it requires a QgsPointXY() as input:

nearestNeighbor(self, point: QgsPointXY, neighbors: int = 1, maxDistance: float = 0)

which makes me think it does not support 3D distances. However, you can also hand over a QgsGeometry() instead:

nearestNeighbor(self, geometry: QgsGeometry, neighbors: int = 1, maxDistance: float = 0) -> List[int] Returns nearest neighbors to a geometry. The number of neighbors returned is specified by the neighbors argument.

QgsGeometry() does support Z values. However, for QgsGeometry() there is (as far as I know) no method to measure 3D distances. For example QgsGeometry().distance() says:

QgsGeometry objects are inherently Cartesian/planar geometries, and the distance returned by this method is calculated using strictly Cartesian mathematics.

So for example

geom1 = QgsGeometry.fromWkt('Point(1 1 0')
geom2 = QgsGeometry.fromWkt('Point(1 1 5')
print(geom1.distance(geom2))

returns 0.0.

Additionally, you cannot use a QgsPoint().

There is no more information about it. So taking a look at QgsSpatialIndexKDBush() it says:

A very fast static spatial index for 2D points based on a flat KD-tree.

Which makes me think "why is 2D mentioned here? Is it because only the KDBush index does not support 3D?"

Is there any documentation about 3D support or a way to figure out if and how QgsSpatialIndex() measures (3D) distances? Or an alternative index supporting 3D measures?

Answer 1

I could not find any documation on this, but according to my tests, no, QgsSpatialIndex() does not support 3D measures. Try:

vl = QgsVectorLayer("PointZ", "temp", "memory") # adding and saving to a layer, otherwise feature ids are all the same
pr = vl.dataProvider()
vl.startEditing()
f1 = QgsFeature()
geom1 = QgsGeometry.fromWkt('Point(1 1 0')
f1.setGeometry(geom1)
pr.addFeature(f1)
f2 = QgsFeature()
geom2 = QgsGeometry.fromWkt('Point(1 1 10')
f2.setGeometry(geom2)
pr.addFeature(f2)
f3 = QgsFeature()
geom3 = QgsGeometry.fromWkt('Point(1 1 500')
f3.setGeometry(geom3)
pr.addFeature(f3)
f4 = QgsFeature()
geom4 = QgsGeometry.fromWkt('Point(1 2 0')
f4.setGeometry(geom4)
pr.addFeature(f4)
vl.commitChanges()
vl.updateExtents() 
QgsProject.instance().addMapLayer(vl)
idx = QgsSpatialIndex(flags=QgsSpatialIndex.FlagStoreFeatureGeometries)
#idx = QgsSpatialIndex() # there is no difference with or without stored geometries
idx.addFeature(f1)
idx.addFeature(f2)
idx.addFeature(f3)
idx.addFeature(f4)
nn = idx.nearestNeighbor(geom1,-1)
print(nn)
for fid in nn:
    print('Feat ' + str(fid) + ' WKT: ' + idx.geometry(fid).asWkt())

returns:

[1, 3, 2, 4]
Feat 1 WKT: PointZ (1 1 0)
Feat 3 WKT: PointZ (1 1 500)
Feat 2 WKT: PointZ (1 1 10)
Feat 4 WKT: PointZ (1 2 0)

while it actually should return:

[1, 4, 2, 3]
Feat 1 WKT: PointZ (1 1 0)
Feat 4 WKT: PointZ (1 2 0)
Feat 2 WKT: PointZ (1 1 10)
Feat 3 WKT: PointZ (1 1 500)

Please correct me if my tests are invalid, incorrect, I am missing something or some docs state something different.