Iterating over two very large dataframes: finding the closest LineString (1/150k) to a Point (1/1.2m)

I have a dataset of calls (1.2m rows), each with a 'Point' geometry. Likewise I have a street dataset (150k rows), each with a LineString geometry. My goal is to identify, for each call, the closest street.

I have the following implementation (from Tenkanen, Heikinheimo & Aagesen), which works very efficiently, however relies on the prior calculation of each LineString centroid. i.e. I'm finding the closest centroid (Point) to each call Point.

I would prefer to not have to do make this assumption and instead actually find the closest street LineString to each call Point. Does anyone know how I could modify the code below to achieve such an objective?

def get_nearest(src_points, candidates, k_neighbors=1):
"""Find nearest neighbors for all source points from a set of candidate points"""

# Create tree from the candidate points
tree = BallTree(candidates, leaf_size=15, metric='haversine')

# Find closest points and distances
distances, indices = tree.query(src_points, k=k_neighbors)

# Transpose to get distances and indices into arrays
distances = distances.transpose()
indices = indices.transpose()

# Get closest indices and distances (i.e. array at index 0)
# note: for the second closest points, you would take index 1, etc.
closest = indices
closest_dist = distances

# Return indices and distances
return (closest, closest_dist)

def nearest_neighbor(left_gdf, right_gdf, return_vals=False):
"""
For each point in left_gdf, find closest point in right GeoDataFrame and return them.
"""

left_geom_col = left_gdf.geometry.name
right_geom_col = right_gdf.geometry.name

# Ensure that index in right gdf is formed of sequential numbers
right = right_gdf.copy().reset_index(drop=True)

# Parse coordinates from points and insert them into a numpy array as RADIANS
left_radians = np.array(left_gdf[left_geom_col].apply(lambda geom: (geom.x * np.pi / 180, geom.y * np.pi / 180)).to_list())
right_radians = np.array(right[right_geom_col].apply(lambda geom: (geom.x * np.pi / 180, geom.y * np.pi / 180)).to_list())

# Find the nearest points
# -----------------------
# closest ==> index in right_gdf that corresponds to the closest point
# dist ==> distance between the nearest neighbors (in meters)

# Return points from right GeoDataFrame that are closest to points in left GeoDataFrame
closest_points = left_gdf # right.loc[closest]

# Ensure that the index corresponds the one in left_gdf
closest_points = closest_points.reset_index(drop=True)

# Add the head and tail node IDs of the closest street
if return_vals:
closest_points['u'] = right.loc[closest,'u'].reset_index(drop=True)
closest_points['v'] = right.loc[closest,'v'].reset_index(drop=True)

return closest_points
• I wonder if the SpatialIndex.nearest method in GeoPandas would help you. In general, you want to query an RTree of the LineStrings to find ones nearby, then filter on distance to find the absolute closest. The SpatialIndex is GoePandas wrapper around some available RTrees. But you could also use shapely.strtree.STRtree, rtree, or whatever pygeos has for RTrees directly. Nov 17 '21 at 1:09