pythonshapely

Shapely unable to split line on point due to precision issues


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.


Solution

  • 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