Shapely ST_Split

split.py

from shapely.geometry import Point, LineString, Polygon, GeometryCollection, MultiLineString
from shapely.ops import polygonize, unary_union, nearest_points
from shapely.wkt import loads as load_wkt

def split_polygon_line(polygon, line):
    """Split a Polygon with a LineString"""
    assert(isinstance(polygon, Polygon))
    assert(isinstance(line, LineString))
    boundary = polygon.boundary
    union = boundary.union(line)
    collection = []
    for geometry in list(polygonize(union)):
        if polygon.contains(geometry.representative_point()):
            collection.append(geometry)
    collection = GeometryCollection(collection)
    return collection

def split_line_line(line1, line2):
    """Split a LineString with another LineString"""
    if isinstance(line2, Polygon):
        line2 = line2.boundary
    assert(isinstance(line, LineString))
    assert(isinstance(line2, LineString))
    assert(not line.relate_pattern(line2, '1********'))
    return GeometryCollection([l for l in line1.difference(line2)])

def split_line_point(line, point):
    p1, p2 = nearest_points(line, point)
    assert(p1.distance(p2) == 0.0)
    coords = line.coords
    for n in range(0, len(coords)-1):
        # http://stackoverflow.com/a/328122/1300519
        a = p1
        b = coords[n]
        c = coords[n+1]
        crossproduct = (c.y - a.y) * (b.x - a.x) - (c.x - a.x) * (b.y - a.y)
        if abs(crossproduct) > epsilon :
            continue
        dotproduct = (c.x - a.x) * (b.x - a.x) + (c.y - a.y)*(b.y - a.y)
        if dotproduct < 0 :
            continue
        squaredlengthba = (b.x - a.x)*(b.x - a.x) + (b.y - a.y)*(b.y - a.y)
        if dotproduct > squaredlengthba:
            continue
    print(n)

polygon = Polygon([(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)])
line = LineString([(5, -20), (5, 20)])

collection = split_polygon_line(polygon, line)

union = unary_union(collection)

assert(polygon.area == collection.area)
assert(polygon.area == union.area)
assert(polygon.symmetric_difference(union).empty)

line2 = polygon.boundary

collection = split_line_line(line, line2)

assert(collection == load_wkt('GEOMETRYCOLLECTION (LINESTRING (5 -20, 5 0), LINESTRING (5 0, 5 10), LINESTRING (5 10, 5 20))'))

collection = split_line_point(line, Point(5, 0))

print(collection)

Dumping this here so I don't lose it, but it should eventually make it's way into Shapely.ops.