Bulk loading into pyrtree index

Question

I'm using pyrtree to create an index for a point in polygon operation. I am currently trying to insert 1.8 million rows into the index and it takes more than an hour to do so. I am unable to use rtree in pyspark and hence chose to use pyrtree. I read about bulk inserting into rtree but I did not find anything similar in pyrtree. Below is the script I am using for the insert (single insert)

from shapely.geometry import Point, shape
from pyproj import Proj, transform
from pyrtree import RTree,Rect
places_idx_new=RTree()
for index,row in places_data.iterrows():
    pt=Point(transform(Proj(init='EPSG:4326'),Proj(init='EPSG:3857'),row['longitude'],row['latitude']))
    pt_buff=pt.buffer(places_data['category_buffer'][index])
    places_idx_new.insert(index,Rect(*pt_buff.bounds))

Please let me know if there is a way to bulk load into pyrtree or if there is a faster way to do this and If I can provide any further information on this

Answer 1

This is making a big assumption: What you're trying to do is find all the "places" that are within their individual category_buffer distance of a new point (or new points).

A second assumption is that you have access to scipy on the cluster you're working with for pyspark (though I think that's a safer call to make).

At any rate assuming this is correct you can do it by extrapolating a third dimension. This third dimension will still allow the distances on the 2d plane to be the same, but will set the height of each "place" so that the distance in the 3d dataset is the same for each point such that a right angle triangle with a base the length of the category_buffer (in the same dimension as the original 2d plane) will have a hypotenuse the length of the longest buffer and the height of the triangle will be used as the third dimension. e.g.

In the image above if you search Point 1 and Point 2 for points within 5 units in the 2d plane (x axis in the image) then you will get both Place A and Place B, when really you're expecting only Place B. If you add the height factor and search the same distance, Point 1 will be within buffer distance of Place B, and only Point 2 will be within buffer distance of Place A and B.

Why bother? If you can use a single search distance this kind of search is simple to implement using a KDTree. In this case we'll use scipy.spatial.cKDTree to do the job. I'm going to assume you have 3 numpy arrays: one of eastings, one of northings (pyproj can project entire arrays in one hit rather than individual points if it helps) and one of the category_buffer.

import numpy as np

max_buffer = np.max(category_buffer)
z = np.sqrt((max_buffer ** 2) - (category_buffer ** 2))
xyz = np.column_stack([eastings, northings, z])

import scipy.spatial

tree = scipy.spatial.cKDTree(xyz)

At this point you can query point at a time (by setting the 3rd dimension to 0), or by an array of points (again, with the 3rd dimension at 0).

idx = tree.query_ball_point([point.x, point.y, 0], max_buffer)

This returns a list of indexes of places in range, which you should be able to match up to your original places_data pretty easily. Querying multiple points will get you a list of lists.