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.

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?

share|improve this question

1 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.

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.