Saturday 19 December 2015

rtree python polygon index


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?




No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...