1

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?

1 Answer 1

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.

1
  • Same goes for overlay_nearest() function; you can test with array_to_string(overlay_nearest(@layer,geom_to_wkt($geometry),limit:=-1))
    – MrXsquared
    Jan 24 at 20:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.