snap points to nearest OSMnx edges.

snap_point_to_osmnx_edge.py

# snap points to nearest OSMnx edges.
import argparse
import geopandas as gpd
import osmnx as ox
import pandas as pd

parser = argparse.ArgumentParser()
parser.add_argument("point_file", help="point list csv file")
parser.add_argument("-x", "--longitude", help="column name for longitude", default="longitude")
parser.add_argument("-y", "--latitude", help="column name for latitude", default="latitude")
parser.add_argument("-b", "--bbox", help="comma separated bounding box for osmnx graph_from_bbox(). (default: total_bounds of point_file)")
parser.add_argument("-n", "--network", help="network type for osmnx graph_from_bbox()", default="drive")
parser.add_argument("-o", "--output_file", help="output filename", default="points_on_edge.csv")
args = parser.parse_args()

df = pd.read_csv(args.point_file, encoding="utf-8")
points_gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df[args.longitude], df[args.latitude]), crs=ox.settings.default_crs)

if args.bbox:
    bbox = [float(x) for x in args.bbox.split(",")]
else:
    bbox = points_gdf.total_bounds
print("bbox", bbox)
G = ox.graph_from_bbox(bbox, network_type=args.network, simplify=False)
Gp = ox.project_graph(G)
points = points_gdf.to_crs(Gp.graph["crs"])
ne = ox.distance.nearest_edges(Gp, points.geometry.x, points.geometry.y)

edges_gdf = ox.graph_to_gdfs(Gp, nodes=False)
nearest_edges_gdf = edges_gdf.loc[ne].reset_index(drop=True)

# nearest point on edge
# https://spaceandtim.es/code/closest_edges/
project2edge = nearest_edges_gdf.project(points, normalized=True)
points_on_edge = nearest_edges_gdf.interpolate(project2edge, normalized=True)

pointslatlon = gpd.GeoDataFrame(geometry=points_on_edge.to_crs(ox.settings.default_crs))
points_gdf[args.longitude] = pointslatlon.geometry.x
points_gdf[args.latitude] = pointslatlon.geometry.y

points_gdf.drop(columns=["geometry"]).to_csv(args.output_file, index=False)