3

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?

enter image description here

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 enter image description here

Below you can see two of the 16 line obtained splitting line_orange with the polygon: enter image description here enter image description here

1 Answer 1

3

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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