Shapely strange splits when splitting LineString and Polygon

Question

I noticed a strange behavior while splitting a LineString with e Polygon: in the below image I have a Polygon (in green) and two lines. If I split the blue Linestring with the polygon the result of the split is a list of 3 LineString (as expected), but if I split the orange Linestring with the polygon the result is a list of 16 LineString. Instead, I expected only 4 linestring; what am I missing?

The polygon is the result of line_blu.buffer(<value>).intersect(line_orange).buffer(<value>) so I am sure that all the orange point you see in the green area are actually in that polygon.
Just as prove of that if I calculate the points outside the polygon with the following code

xy = [xy for xy in line_orange.coords if not polygon.contains(Point(xy))]

The result is

Below you can see two of the 16 line obtained splitting line_orange with the polygon:

linestring polygon python shapely splitting

Jeremiah England · Accepted Answer

It looks like the orange linestring crosses itself, which is related to issues #600 and #1068 in the shapely github repository. The problem is that shapely's split function splits complex (i.e. self-intersecting) linestrings at their self-intersection points in addition to the points where they intersect the splitter geometry.
Here is a solution that I cooked up for one of those issues which seems to work. Let me know if you have any issues with it or improvements!
""""
Split a complex linestring using shapely.

Inspired by https://github.com/Toblerity/Shapely/issues/1068
"""
from shapely.geometry import LineString, GeometryCollection, Point, Polygon
from shapely.ops import split, snap

def complex_split(geom: LineString, splitter):
    """Split a complex linestring by another geometry without splitting at
    self-intersection points.

Parameters
    ----------
    geom : LineString
        An optionally complex LineString.
    splitter : Geometry
        A geometry to split by.

Warnings
    --------
    A known vulnerability is where the splitter intersects the complex
    linestring at one of the self-intersecting points of the linestring.
    In this case, only the first path through the self-intersection will
    be split.

Examples
    --------
    >>> complex_line_string = LineString([(0, 0), (1, 1), (1, 0), (0, 1)])
    >>> splitter = LineString([(0, 0.5), (0.5, 1)])
    >>> complex_split(complex_line_string, splitter).wkt
    'GEOMETRYCOLLECTION (LINESTRING (0 0, 1 1, 1 0, 0.25 0.75), LINESTRING (0.25 0.75, 0 1))'

Return
    ------
    GeometryCollection
        A collection of the geometries resulting from the split.
    """
    if geom.is_simple:
        return split(geom, splitter)
    
    if isinstance(splitter, Polygon):
        splitter = splitter.exterior

# Ensure that intersection exists and is zero dimensional.
    relate_str = geom.relate(splitter)
    if relate_str[0] == '1':
        raise ValueError('Cannot split LineString by a geometry which intersects a '
                         'continuous portion of the LineString.')
    if not (relate_str[0] == '0' or relate_str[1] == '0'):
        return geom

intersection_points = geom.intersection(splitter)
    # This only inserts the point at the first pass of a self-intersection if
    # the point falls on a self-intersection.
    snapped_geom = snap(geom, intersection_points, tolerance=1.0e-12)  # may want to make tolerance a parameter.
    # A solution to the warning in the docstring is to roll your own split method here.
    # The current one in shapely returns early when a point is found to be part of a segment.
    # But if the point was at a self-intersection it could be part of multiple segments.
    return split(snapped_geom, intersection_points)

if __name__ == '__main__':
    complex_line_string = LineString([(0, 0), (1, 1), (1, 0), (0, 1)])
    splitter = LineString([(0, 0.5), (0.5, 1)])

out = complex_split(complex_line_string, splitter)
    print(out)
    assert len(out) == 2

# test inserting and splitting at self-intersection
    pt = Point(0.5, 0.5)
    print(f'snap: {snap(complex_line_string, pt, tolerance=1.0e-12)}')
    print(f'split: {split(snap(complex_line_string, pt, tolerance=1.0e-12), pt)}')

Note: I have not tested with Polygons, only LineStrings. But I just added these two lines which should make it work for Polygons.
    if isinstance(splitter, Polygon):
        splitter = splitter.exterior

Shapely strange splits when splitting LineString and Polygon

One Answer

Add your own answers!

Ask a Question