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:

Answer 1

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