How can I use R-Tree to find points within a distance in Spatialite?

Question

Create a table with 3D points:

create table points (coords pointz);
select CreateSpatialIndex('points', 'coords');

Insert 13k rows:

insert into points values (geometryfromtext('pointz(...)'))

Do

SELECT
    p.rowid, count(*)
FROM
    points p, points v
WHERE
    distance(v.coords, p.coords) < 1
GROUP BY
    p.rowid;

Will see an output of several (1..10) points per second (will take a few hours to process 13K points). EXPLAIN shows that the index is not used. A query with PTDistanceWithin function works the same way.

How can I use the spatial index? I checked if it's possible to use an MBR, but those functions are for 2D only.

added: here's a function reference, and all the functions I could use in usual case, work in 2D only!

Answer 1

Here are two non-Spatialite solutions.

For Rtree with Python, try this example with 13000 points -- should take a few seconds:

from random import randrange
from rtree import index
from math import sqrt

# Create a 3D index
p = index.Property()
p.dimension = 3
idx3d = index.Index(properties=p)

# Make and index random data
coords = []
for id in range(13000):
    coord = (randrange( 100000,  500000),
             randrange(1000000, 5000000),
             randrange(      0,    5000))
    coords.append(coord)
    idx3d.add(id, coord)

# Find closest pair for the first 10 points
for id1 in range(10):
    nearest = list(idx3d.nearest(coords[id1], 2))
    assert id1 == nearest[0]
    id2 = nearest[1]
    c1 = coords[id1]
    c2 = coords[id2]
    # Pythagorean theorem
    dist = sqrt(sum([(a - b)**2 for a, b in zip(c1, c2)]))
    print '%i <-> %i : %.1f'%(id1, id2, dist)

For PostGIS 2.0/PostgreSQL 9.1, there is a new KNN operator <-> to do a similar nearest search (but I'm not sure if it is 3D aware). Or there is ST_3DDwithin to get points in a 3D radius search. See these answers to get a further idea:

• http://stackoverflow.com/a/7947019/327026
• http://gis.stackexchange.com/a/21905/1872

Answer 2

You should use the 2D version and the RTree index. Use of index is not automatic - you have to explicitly invoke (i.e. include it in the SQL query)

The SpatiaLite cookbook explains the background to the RTree index at http://www.gaia-gis.it/gaia-sins/spatialite-cookbook/html/rtree.html.

There are a lot of examples in the Cookbook, especially in http://www.gaia-gis.it/gaia-sins/spatialite-cookbook/html/neighbours.html.