Search code examples
pythongeometryshapely

Split a linestring at each point it crosses a polygon in shapely


I have a linestring which spans various polygons, stored as GeoJsons. I want to split the line into individual parts within each polygon area. I have not been able to achieve this yet, however. This is a reproducible example for what I have so far:

from shapely.geometry import Point, LineString, Polygon
from shapely.ops import split
from matplotlib import pyplot as plt

poly = {
    'type': "Feature",
    'geometry': {
        "type": "Polygon",
        "coordinates": ([(2,2),(2,4),(4,4),(4,2)]),
    },
    'properties': {
        'id': 'A'
    }
}

line = {
    'type': "Feature",
    'geometry': {
        "type": "Linestring",
        "coordinates": ([(3,3),(5,1)]),
    },
    'properties': {
        'id': 'A'
    }
}

poly = Polygon(poly['geometry']['coordinates'])
line = LineString(line['geometry']['coordinates'])

x,y = poly.exterior.xy
x2,y2 = line.coords

plt.plot(x, y, x2, y2)
plt.show()

This code produces the following square polygon with a linestring crossing it: enter image description here

I've then tried to split the line via the polygon, as so:

new_lines = split(line, poly)

But I get the following output which doesn't seem correct:

GEOMETRYCOLLECTION (LINESTRING (3 3, 4 2), LINESTRING (4 2, 5 1))

I was expecting three lines, one existing inside the square polygon, and then two existing separately outside the polygon.


Solution

  • It looks like split is working as you expect, but the line's coordinates are plotted in a misleading way - should use .xy to get the line xs and ys as with the polygon:

    from shapely.geometry import Point, LineString, Polygon
    from shapely.ops import split
    from matplotlib import pyplot as plt
    
    poly = Polygon([(2,2),(2,4),(4,4),(4,2)])
    original_line = LineString([(3,3),(5,1)])
    line = LineString([(3,1), (3,5)])
    x, y = poly.exterior.xy
    x2, y2 = line.xy
    x3, y3 = original_line.xy
    
    plt.plot(x, y, x2, y2)
    plt.show()
    
    plt.plot(x, y, x3, y3)
    plt.show()
    

    image with original coords

    image with expected coords

    >>> print(split(line, poly))
    GEOMETRYCOLLECTION (LINESTRING (3 1, 3 2), LINESTRING (3 2, 3 4), LINESTRING (3 4, 3 5))