1

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

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

0

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.