Is it possible to fit a coordinate to a street in OSMnx?

Question

OSMnx provides solution to calculate the shortest path between two nodes, but I would like to the same with points on streets (I have GPS coordinates recorded from vehicles). I know there is also a method to get the closest node, but I have two question for this problem of mine.

i) When closest node computed is the street where the point is also taken into consideration? (I assume not) ii) If I wanted to implement something like this, I like to know how a street (edge) is represented as a curve (Bézier curve maybe?). Is it possible to get the curve (or the equation of the curve) of an edge?

I asked this question here, because the guidelines for contributing of OSMnx asked it.

Answer 1

OSMnx was recently updated since there have been a couple of requests in this direction (see https://github.com/gboeing/osmnx/pull/234 and references therein). So in the last update, you'll find a function like this:

ox.get_nearest_edge(G, (lat, lon))

It will give you the ID of the nearest edge, which is much better than nearest nodes. However, I think it is more useful to also get the actual distance of the nearest edge in order to check whether or not your data point is on the road or a few thousand meters apart...

To do this, I followed the implementation from https://github.com/gboeing/osmnx/pull/231/files

# Convert Graph to graph data frame
gdf = ox.graph_to_gdfs(G, nodes=False, fill_edge_geometry=True)
# extract roads and some properties
roads = gdf[["geometry", "u", "v","ref","name","highway","lanes"]].values.tolist()
# calculate and attach distance
roads_with_distances = [(road, ox.Point(tuple(reversed((lat,lon)))).distance(road[0])) for road in roads]

# sort by distance
roads_with_distances = sorted(roads_with_distances, key=lambda x: x[1])
# Select closest road
closest_road = roads_with_distances[0]
# Check whether you are actually "on" the road
if closest_road[1] < 0.0001: print('Hit the road, Jack!')

I have the impression that a distance on the order of $10^{-5}$ means that the coordinate is actually "on" the road.

Answer 2

Streets and node in OSMnx are shapely.geometry.LineString, and shapely.geometry.Point objects, so there is no curve, only sequence of coordinates. The technical term for what you described is Map Matching. There are different ways of map matching, the simplest one being geometric map matching in which you find nearest geometry (node or edge) to the GPS point. point to point map matching can be easily achieved using built-in osmnx function ox.get_nearest_node(). If you have a luxury of dense GPS tracks, this approach could work reasonably good. For point to line map matching you have to use shapely functions. The problem with this approach is that it is very slow. you can speed up the algorithm using spatial index, but still, it will not be fast enough for most purposes. Note that geometric map matching are least accurate among all approaches. I wrote a function a few weeks ago that does simple point to line map matching using edge GeoDataFrame and node GeoDataFrame that you can get from OSMnx. I abandoned this idea and now I am working on a new algorithm (hopefully much faster), which I will publish on GitHub upon completion. Meanwhile, this may be of some help for you or someone else, so I post it here. This is an early version of abandoned code, not tested enough and not optimized. give it a try and let me know if it works for you.

def GeoMM(traj, gdfn, gdfe):
"""
performs map matching on a given sequence of points

Parameters
----------

Returns
-------
list of tuples each containing timestamp, projected point to the line, the edge to which GPS point has been projected, the geometry of the edge))

"""

traj = pd.DataFrame(traj, columns=['timestamp', 'xy'])
traj['geom'] = traj.apply(lambda row: Point(row.xy), axis=1)
traj = gpd.GeoDataFrame(traj, geometry=traj['geom'], crs=EPSG3740)
traj.drop('geom', axis=1, inplace=True)

n_sindex = gdfn.sindex

res = []
for gps in traj.itertuples():
    tm = gps[1]
    p = gps[3]
    circle = p.buffer(150)
    possible_matches_index = list(n_sindex.intersection(circle.bounds))
    possible_matches = gdfn.iloc[possible_matches_index]
    precise_matches = possible_matches[possible_matches.intersects(circle)]
    candidate_nodes = list(precise_matches.index)

    candidate_edges = []
    for nid in candidate_nodes:
        candidate_edges.append(G.in_edges(nid))
        candidate_edges.append(G.out_edges(nid))

    candidate_edges = [item for sublist in candidate_edges for item in sublist]
    dist = []
    for edge in candidate_edges:
        # get the geometry
        ls = gdfe[(gdfe.u == edge[0]) & (gdfe.v == edge[1])].geometry
        dist.append([ls.distance(p), edge, ls])

    dist.sort()
    true_edge = dist[0][1]
    true_edge_geom = dist[0][2].item()
    pp = true_edge_geom.interpolate(true_edge_geom.project(p)) # projected point
    res.append((tm, pp, true_edge, true_edge_geom))


    return res