Finding nearest OSM link (road) to GPS points in Python

Question

I have a large dataset (about 50,000,000 rows) of GPS x,y points related to GPS equipped taxis with timestamps and speeds. I also have a shp. of the OSM links and I want to find the nearest link to each GPS point in order to calculate average and stdev of speeds for each link at 5-min intervals.

I am using the following code in Spyder to check the results for the first row of my GPS x,y point dataset and it returns the result (OSM link ID) correctly.

subset_point = points_gdf.head(1).copy()
subset_point.crs = 'EPSG:4326'
subset_point['geometry'] = gpd.points_from_xy(subset_point['lon'], subset_point['lat'])
def find_nearest_link(point, links, search_radius=0.00001):
    nearest_link = links.loc[links.geometry.apply(lambda geom: point.distance(geom) < search_radius)]
    if not nearest_link.empty:
        return nearest_link.iloc[0]['osm_id']
    else:
        return None
subset_point['nearest_Link'] = subset_point['geometry'].apply(lambda point: find_nearest_link(point, osm_links_3))
print(subset_point[['lon', 'lat', 'nearest_Link']])

However, when performing the following code for all my GPS x,y points the procedure becomes very slow:

def find_nearest_link(point, links, search_radius=0.00001):
    nearest_link = links.loc[links.geometry.apply(lambda geom: point.distance(geom) < search_radius)]
    if not nearest_link.empty:
        return nearest_link.iloc[0]['osm_id']
    else:
        return None
points_gdf['nearest_Link'] = points_gdf['geometry'].progress_apply(lambda point: find_nearest_link(point, osm_links_3))
print(points_gdf[['lon', 'lat', 'nearest_Link']]).

Should I use another method? Is there something more appropriate given the length of my dataset?

Answer 1

The problem with your code is that it scales badly with the number of points, and for each point it will check the distance to every single links (links.geometry.apply(lambda geom: point.distance(geom) < search_radius)).

A simple solution with indexing will result in a massive speedup. Using geopandas.sjoin_nearest() will process your calcul extremely quickly:

import geopandas as gpd
points_gdf = gpd.read_file("points.gpkg", layer='points')
links_gdf = gpd.read_file("links.gpkg", layer='links')
points_joined = gpd.sjoin_nearest(points_gdf, links_gdf, distance_col="dist", how='left')

Profiling this solution versus yours gives the following for a link layer of 1,000,000 entities:

• for 5 points:
  • find_nearest_link: 28.7s (5.75s by point)
  • sjoin_nearest: 234.7 ms = 0.235s
• for 20 points:
  • find_nearest_link: 1m53s (5.68s by point)
  • sjoin_nearest: 236.6ms = 0.237s
• for 100 points:
  • find_nearest_link: 9m59s (6.0s by point)
  • sjoin_nearest: 237.0ms = 0.237s
• for 50,000 points:
  • find_nearest_link: NOT TESTED but mathematically (considering 5.5s by point like for the small samples / very conservative): ~76h
  • sjoin_nearest: 1.787s

As you can see for sjoin_nearest(), the amount of time taken slowly increases with the number of points used. With only an x7 increase in time for an x5,000 increase in the number of points.

While your function scales linearly with the amount of points.

This shows the difference in scalability that can exist for algorithms that are built with a linear time complexity (O(n)) vs logarithm time complexity (O(log(n))

PS: I don't know how many links you have in your OSM Shapefile, but it will still probably take a few minutes for 50,000,000 points.