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.
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