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

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!

• Are you willing to use PostGIS? Or use the Rtree python package, (i.e., without a database)? Both PostGIS and Rtree support 3D R-tree indices without hassle and are very fast. – Mike T Mar 21 '12 at 21:30
• Can you, please, write an example, how to find all pairs of points close one to another? I can use it too. I'd prefer PostGIS over SpatiaLite after learning the latter more. After installing Python library for spatialite, compiling GEOS, P-something, and seeing things working differently in Python than in spatialite shell, I want SpatiaLite go to hell. – culebrón Mar 21 '12 at 21:41

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)

# Find closest pair for the first 10 points
for id1 in range(10):
nearest = list(idx3d.nearest(coords[id1], 2))
assert id1 == nearest
id2 = nearest
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:

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.