Create a table with 3D points:

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

Insert 13k rows:

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


    p.rowid, count(*)
    points p, points v
    distance(v.coords, p.coords) < 1

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
  • 1
    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))
    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:


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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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