Snap points shapefile to line shapefile using shapely

Question

I have 2 datasets: 1) a road network shapefile which has many polylines and 2) a bus stop shapefile which has many points representing bus stops. My points shapefile has CRS of EPSG 4326 and my line shapefile has CRS of EPSG 29617. I am not sure if I should and how to convert them.

I would like to use Python only and available packages (e.g. shapely, fiona, ogr) to snap the closest bus stop that is within 45 metres to the nearest road. Please note and see picture: there are multiple bus stops (green dots) near a road(brown line). Only the closest green dot (i.e. one green dot only) must be snapped to the brown line.

I have written this code based on many forums which have parts of the answer. I have never coded in Python before.

# path of inputs 
line = "C:\folder\line.shp"
point = "C:\folder\point.shp"

# open line and points data
gdf_segments = geopandas.read_file(line)
gdf_hec_points = geopandas.read_file(point)

#join all line and points
shply_line = gdf_segments.geometry.unary_union
shply_point = gdf_hec_points.geometry.unary_union

# perform interpolation and project point to line
pt_interpolate = shply_line.interpolate(shply_line.project(shply_point))

The result of the code should be a new point shapefile where the points have been moved to the line.

Reference: I found the tool in ArcGIS (http://pro.arcgis.com/en/pro-app/tool-reference/editing/snap.htm) but I cannot use any software like QGIS, GRASS, ArcGIS etc. The code needs to be in Python only.

Answer 1

The methodology is called linear referencing and a solution was given by Mike T in Coordinate of the closest point on a line with Shapely. There is also a recipe in the Python Geospatial Analysis Cookbook (Snapping a point to the nearest line )

"This super common spatial task is for all the GPS junkies who want their GPS coordinates to snap to an existing road" ...

The problem:

import geopandas as gpd
gdf_segments = gpd.read_file("line.shp")
shply_line = gdf_segments.geometry.unary_union
point  = gpd.read_file('points.shp')
point.crs
{'init': 'epsg:4326'}
# reproject the points
point = point.to_crs(gdf_segments.crs)
print(point)
   id                                      geometry
0   1  POINT (165.2232667307835 -581.6023098314181)
1   2  POINT (458.0395231332805 -626.3180927932262)
2   3  POINT (807.1111194563855 -509.4800791162997)
3   4  POINT (1019.150477463845 -1181.659268662333)
4   5  POINT (74.34925616221153 -244.0702704396099)
5   6  POINT (19.53636085704784 -383.9873978914693)

Solution with linear referencing:

for i in range(len(point)):
    print(shply_line.interpolate(shply_line.project( point.geometry[i])).wkt)

POINT (158.2568242503091 -613.1561963606259)
POINT (433.554325720973 -616.258892163531)
POINT (828.8651101528822 -533.8813967264118)
POINT (981.8579545397143 -1193.379775867061)
POINT (74.34925616221147 -233.4585492227979)
POINT (18.88030104343747 -405.9654016474183)

New GeoDataFrame with results:

result = point.copy()
result['geometry'] = result.apply(lambda row: shply_line.interpolate(shply_line.project( row.geometry)), axis=1)
print(result)
   id                        geometry
0   1  POINT (158.2568242503091 -613.1561963606259)
1   2  POINT (433.554325720973 -616.258892163531)
2   3  POINT (828.8651101528822 -533.8813967264118)
3   4  POINT (981.8579545397143 -1193.379775867061)
4   5  POINT (74.34925616221147 -233.4585492227979)
5   6  POINT (18.88030104343747 -405.9654016474183)
result.to_file("new_points.shp")

NEW

Buffer the lines (45m)

buff = shply_line.buffer(45)

Select the points within 45m (point in polygon):

from geopandas.tools import sjoin
buff = gpd.GeoDataFrame( geometry=[buff])
pointInPolys = sjoin(point,buff, how='left')
pointInPolys
   id                  geometry                      index_right
0   1  POINT (165.2232667307835 -581.6023098314181)          NaN
1   2  POINT (458.0395231332805 -626.3180927932262)          0.0
2   3  POINT (807.1111194563855 -509.4800791162997)          NaN
3   4  POINT (1019.150477463845 -1181.659268662333)          NaN
4   5  POINT (74.34925616221153 -244.0702704396099)          0.0
5   6  POINT (19.53636085704784 -383.9873978914693)          0.0
point45 = pointInPolys.dropna()
point45
   id                  geometry                      index_right
1   2  POINT (458.0395231332805 -626.3180927932262)          0.0
4   5  POINT (74.34925616221153 -244.0702704396099)          0.0
5   6  POINT (19.53636085704784 -383.9873978914693)          0.0
result2 =  point45.copy()
result2['geometry'] = result2.apply(lambda row:    shply_line.interpolate(shply_line.project( row.geometry)), axis=1)
print(result2)
   id                  geometry                      index_right
1   2  POINT (433.554325720973 -616.258892163531)          0.0
4   5  POINT (74.34925616221147 -233.4585492227979)          0.0
5   6  POINT (18.88030104343747 -405.9654016474183)          0.0

Points in Green