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[0]
    closest_dist = distances[0]

    # 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)

    closest, dist = get_nearest(src_points=left_radians, candidates=right_radians)

    # 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.