shapely wrong output -- when snapping points to line shapefile

Question

I have a series of nodes to be connected to a road network, so as to calculate the cost of transportation. I used shapely to snap the nodes to the road network polylines. Most of the connections worked well, however, there are several connecting to very far road segments, instead of a near road. I don't understand what could have gone wrong here.

strange snapping results

here is the essential part of the codes:

# points are the nodes
# roadLine is the road network in linestring
ptsPrj = points.copy() # find the points project to locations on the road network 
ptsPrj['geometry'] = ptsPrj.apply(lambda row: roadLine.interpolate(roadLine.project(row.geometry)), axis = 1)
# join the points with its projecting points on the road network
joinPts = points.copy()
for i in range(len(points)):
 joinPts.at[i, 'geometry'] = LineString([points.at[i,'geometry'], ptsPrj.at[i, 'geometry']])

Answer 1

For anyone who might be interested in the future, I found a solution:

# points are the nodes; roadLine is the road network in linestring
join = points.copy()
join['geometry'] = join.apply(lambda row: nearest_points(row.geometry, roadLine), axis = 1)
join['line'] = join.apply(lambda row: LineString(row.geometry), axis = 1)
join = join.set_geometry('line')

So, essentially, using 'nearest_points' appears to be more reliable than using the interpolate-project functions to find shortest path.

the corrected output is as the graph below. corrected graph