Shapely floating problems with split()

Based on the answer of sgillies I want to use Shapely's interpolate() to find the points at which to split. But I experience that split() is sensitive to the points being precisely on the geometry you want to split. How to cope with these Floating point issues.

I introduce the following function which should return line segment of equal length:

def get_linesegments(line, n):
points = MultiPoint([line.interpolate(i/n, normalized=True) for i in range(1, n)])
print(points.wkt)
return split(line, points)

Now generate a line:

line = LineString([(0,1), (1,0), (2, 1), (1, 2)])

When the line in split in three parts, the function works correct:

>> lines = get_linesegments(line, 3)
MULTIPOINT (1 0, 2 1)
>> lines.wkt
'GEOMETRYCOLLECTION (LINESTRING (0 1, 1 0), LINESTRING (1 0, 2 1), LINESTRING (2 1, 1 2))'

But when the line is split into 2 parts, we get floating points and the function returns the complete line:

>> lines = get_linesegments(line, 2)
MULTIPOINT (1.5 0.5000000000000001)
>> lines.wkt
'GEOMETRYCOLLECTION (LINESTRING (0 1, 1 0, 2 1, 1 2))'

Another workaround is to sacrifice some precision and use buffer+difference:

def get_linesegments(line, n):
points = MultiPoint([line.interpolate(i/n, normalized=True) for i in range(1, n)])
return line.difference(points.buffer(1e-13))
>>> get_linesegments(line, 2).wkt
'MULTILINESTRING ((0 1, 1 0, 1.499999999999929 0.4999999999999294), (1.500000000000071 0.5000000000000707, 2 1, 1 2))'

>>> get_linesegments(line, 3).wkt
'MULTILINESTRING ((0 1, 0.9999999999999293 7.071566359159048e-14), (1.000000000000071 7.06582312953894e-14, 1.999999999999929 0.9999999999999294), (1.999999999999929 1.000000000000071, 1 2))'

I created a workaround where the distance is checked with a certain limit, and a list is generated with line segments.

def get_linesegments(line, n):
segments = [line]
points = MultiPoint([line.interpolate(i/n, normalized=True) for i in range(1, n)])
for point in points:
lastline = segments[-1]
for ix, (c1, c2) in enumerate(zip(lastline.coords[:-1], lastline.coords[1:])):
if LineString([c1, c2]).distance(point) < 1e-8:
segments[-1] = LineString(lastline.coords[:ix+1] + [point.coords])
if point.coords == c2:
segments.append(LineString(lastline.coords[ix+1:]))
else:
segments.append(LineString([point.coords] + lastline.coords[ix+1:]))
break
return GeometryCollection(segments)

but... it smells a bit