Shapely floating problems with split()

Question

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

Answer 1

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[0]])
                if point.coords[0] == c2:
                    segments.append(LineString(lastline.coords[ix+1:]))
                else:
                    segments.append(LineString([point.coords[0]] + lastline.coords[ix+1:]))
                break
    return GeometryCollection(segments)

but... it smells a bit

Answer 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))'