2

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.

enter image description here

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.

2

This question is duplicated in stackoverflow: snap points shapefile to line shapefile using shapely

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:

enter image description here

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")

enter image description here

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,new, 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['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

enter image description here

  • how do i apply the 45m snapping distance to the interpolate function? Only the closest point within 45m can be snapped to a line segment. Therefore, only 1 point can be snapped to a line segment. – poli science Dec 25 '18 at 21:12
  • Look at New in the answer. – gene Dec 26 '18 at 10:31
  • i edited my question with a picture to show multiple bus stops (green) near a road (brown). Only the closest bus stop within the 45m buffer must be snapped. There should only be 1 green dot on the brown line. – poli science Dec 26 '18 at 16:58
  • I'm sorry, but this can be done by yourself from the solution, It's easy – gene Dec 28 '18 at 10:50
  • This code is very useful for snapping the points. Just a small clarification; I don't see any way to obtain measures of the line file. Is there a way to achieve this? – April 1 hour ago

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.