rtree python polygon index

Question

Sorry for a very basic questions about rtree.

I would like to check whether the point is inside the polygon and the distance (point to the nearest point in polygon) as well. I was thinking about using No-SQL database, but both MongoDb and GeoCouch support Polygon records, but don't provide searches like: - what polygons are around the point (for example, I want to find near what polygons a user is currently located; and after that check what in the precise distance from the user to the nearest polygon point)

So, I want to figure out how to make these kind of spatial searches using python and related libraries only (shapely, geopy, rtree). Here are the steps I am following:

• I have shapefile with the targeted polygons (saved from QGis)
• I want to make Rtree Index, to improve performance when working with Shapely and geopy

I am wondering, whether I can:

• make Rtree Index for shapefile directly?

If not. I can use shapefile lib to collect polygon coordinates or bounding boxes:

>>> shapes[33].points
[[120.9937913, 24.7902669], [120.9939797, 24.790369], [120.994035, 24.7902849], [120.9938467, 24.7901828], [120.9937913, 24.7902669]]

>>> shapes[33].bbox
[120.9937913, 24.7901828, 120.994035, 24.790369]

Now I am a bit confused, and trying to understand how to create rtree index file for my polygons:

>>> idx = index.Index()  //creating index
>>> left, bottom, right, top = (24.7901828, 24.790369, 120.9937913, 120.994035) // can I put bbox coordinates here? I assume that coordinates are in the form [xmin, xmax, ymin, ymax]
>>> idx.insert(0, (left, bottom, right, top)) 

Would it be a right way for me to make Rtree Index?:

import shapefile
from rtree import index 
sf = shapefile.Reader("polygons.shp")
idx = index.Index()
shapes = sf.shapes()
id_num = 0
for shape in shapes:
    sh = shape.bbox
    xmin = sh[1]
    ymin = sh[0]
    xmax = sh[3]
    ymax = sh[2]
    idx.insert(id_num, (xmin,ymin,xmax,ymax))
    id_num = id_num + 1 

After I execute the script above, I can get the following result from idx.bounds:

>>> idx.bounds
[24.7796055, 120.9839184, 24.7986571, 121.0025152]

To what bounds these coordinates are related to? Is it the bbox for the all the polygons?

Is this a right way to check whether there are polygons near the Point? Assuming that the point coordinates are 24.789077, 120.999792

>>> list(idx.nearest((24.789077, 24.789077, 120.999792, 120.999792)))
[1, 2, 13, 29, 61, 62, 59, 60, 27, 28, 55, 5, 26, 54, 78, 49, 50, 3, 9, 10, 22, 72, 43, 39, 42, 85, 40, 8, 17, 76, 77, 36, 73, 71, 33, 34, 69, 63, 66, 64]

I am getting these results and not quite sure what they mean.. What does "nearest" exactly mean here? I cannot measure the distance here, like 5 meters.. so what is the "nearest" output here and how it's calculated?

Answer 1

Your questions are answered in Rtree's docs, which come with the source and are online at http://toblerity.github.com/rtree/tutorial.html.

You can make your shapefile indexing code much simpler:

for i, shape in enumerate(shapes):
    idx.insert(i, shape.bbox)

Use Python's enumerate() function for this kind of thing whenever possible. It's tidy and foolproof.