Shapely unable to split line on point due to precision issues

Question

I am trying to interpolate a point onto a LineString in Shapely and then split the linestring accordingly. However, due to a precision error, Shapely thinks that the interpolated point is not on the linestring, and therefore the split operation does not work.

Here is an example:

from shapely.ops import split
from shapely.geometry import LineString, Point

### Initialize point and line
line = LineString([(0.123,0.456),(5.678,7.890),(12.135,6.789)])    
point = Point(4.785,8.382)   

### Interpolate point onto line
new_point = line.interpolate(line.project(point))
print new_point
>> POINT (5.593949278213755 7.777518800043393)

### BUT: line does not intersect the interpolated point
line.intersects(new_point)
>> False

### EVEN THOUGH: distance between them is essentially (not exactly) zero
line.distance(new_point)
>> 0.0

### THEREFORE: line cannot be split using the new point
len(split(line, new_point))
>> 1

I think the problem is as follows:
1. I rounded the original point/line coordinates so that they do not run up against the precision limits of the machine.
2. However, the INTERPOLATED point has very high precision. I don't know how to control this.
3. In theory, I could round the coordinates of this new point, but that doesn't seem to ensure that the new point is on the line either.

Related questions here, here, and here.

Confirmation

Answer 1

I figured out a somewhat hacky solution. If anyone posts a better one I will accept that instead.

# After the code above: 

### Create a buffer polygon around the interpolated point
buff = new_point.buffer(0.0001)

### Split the line on the buffer
first_seg, buff_seg, last_seg = split(line,buff)

### Stitch together the first segment, the interpolated point, and the last segment 
line = LineString(list(first_seg.coords) + list(new_point.coords) + list(last_seg.coords))

line.intersects(new_point)
>> True

Confirmation

Answer 2

I tried myself with the code above but sometimes it does fail due to a large number of split results...apparently the polygon created breaks up the line in several points, don't know why.

I used instead the snap function, and it seems to work:

    snap(line_to_modify,closest_point,0.01)

The returned value is the modified geometry, with the closest point in it. Firstly, you need to use project and interpolate to find the closest point not in the line. you can veryfy it then with intersect

Bismuthic

Answer 3

i have tried to solve the geopandas "split lines by points" intersection problem with overlay difference and snap.

import geopandas as gpd
from shapely.ops import snap

lines_gdf = gpd.read_file("lines_path")
points_gdf = gpd.read_file("point_path")

points_buf_gdf = points_gdf.copy()
points_buf_gdf.geometry = points_gdf.geometry.buffer(0.0001)
# split
splitted = gpd.overlay(lines_gdf, points_buf_gdf, how='difference', keep_geom_type=True).explode(index_parts=True).reset_index(drop=True)
# fill holes
splitted['geometry'] = snap(splitted['geometry'], points_gdf['geometry'].unary_union, 0.00011)

Abomination