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

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?

• How slow is very slow? because 50,000,000 iterations on geometric operations will take "some" time whatever the way you go about it because if each iteration is 1ms = 0.001s it will still take ~14h. The question is do you need to do it multiple times or can you run it once and save the result to file for future use. And in the latter case do you really care if it takes a day? Commented Nov 29, 2023 at 10:14
• However here it will take a very long time because for every point you check the distance from that point to every line which makes it highly exponential. I would personally look for a better solution. Something using spatial indexing would be my best bet like this or geopandas.sjoin_nearest() Commented Nov 29, 2023 at 11:47

## 1 Answer

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.

• Thank you very much @Louis for your valuable answer. I indeed used the `geopandas.sjoin_nearest()` and also defined a max_distance = 1m and it returned the results (OSM link IDs) in less than half an hour. I was aware of the `geopandas.sjoin()` which comes along with predicate=intersects; contains; within; touches; etc. (but nearest is not an option), but I was not aware of the `geopandas.sjoin_nearest()`. The difference in the amount of time taken is indeed extreme! Commented Nov 30, 2023 at 17:14
• Your welcome, if the answer was satisfactory, please accept the answer Commented Nov 30, 2023 at 21:20