Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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!

share|improve this question
    
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

2 Answers 2

up vote 6 down vote accepted

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:

share|improve this answer

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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