[GIS] Shapely floating problems with split()

linestringpythonshapely

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

Best Answer

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

Related Question