23

I am trying to cut a shapely.geometry.Polygon instance in two parts with two lines. For example, in the code below, polygon is a ring and if we cut it with line1 and line2 we should get two partial rings, one w/ 270 degrees and one with 90 degrees. Would there be a clean way to do this?

from shapely.geometry import Point, LineString, Polygon

polygon = Point(0, 0).buffer(2).difference(Point(0, 0).buffer(1))
line1 = LineString([(0, 0), (3, 3)])
line2 = LineString([(0, 0), (3, -3)])
Georgy
  • 12,464
  • 7
  • 65
  • 73
Yuxiang Wang
  • 8,095
  • 13
  • 64
  • 95

3 Answers3

15

There is a function to split one geometry by another in Shapely as of version 1.6.0 (Aug 2017), so there is no need to roll your own anymore. See the docs for: shapely.ops.split(geom, splitter)

Note that the older answers on this thread were written before the splitting function was in Shapely - they are now effectively obsolete.

Georgy
  • 12,464
  • 7
  • 65
  • 73
Keith Ma
  • 314
  • 3
  • 4
7

Ken Watford answered here about using buffer and difference to do the trick, w/ the drawback of losing a bit of the area. An example code below:

from shapely.geometry import Point, LineString, Polygon

polygon = Point(0, 0).buffer(2).difference(Point(0, 0).buffer(1))
line1 = LineString([(0, 0), (3, 3)])
line2 = LineString([(0, 0), (3, -3)])

line1_pol = line1.buffer(1e-3)
line2_pol = line2.buffer(1e-3)

new_polygon = polygon.difference(line1_pol).difference(line2_pol)

Works for now, and I'd be interested to see whether there is another (potentially w/o losing area) method!

Yuxiang Wang
  • 8,095
  • 13
  • 64
  • 95
  • 2
    I think any exact solution would require to implement the math problem behind from scratch, with significant lines of code. Specially in the general terms you are setting the problem. – eguaio Sep 06 '16 at 11:34
5
from shapely.ops import linemerge, unary_union, polygonize
from shapely.geometry import LineString, Polygon

# Define the Polygon and the cutting line
line = LineString([(-5, -5), (5, 5)])
polygon = Polygon([(-1, -1), (1, -1), (1, 1), (-1, 1)])


def cut_polygon_by_line(polygon, line):
    merged = linemerge([polygon.boundary, line])
    borders = unary_union(merged)
    polygons = polygonize(borders)
    return list(polygons)

def plot(shapely_objects, figure_path='fig.png'):
    from matplotlib import pyplot as plt
    import geopandas as gpd
    boundary = gpd.GeoSeries(shapely_objects)
    boundary.plot(color=['red', 'green', 'blue', 'yellow', 'yellow'])
    plt.savefig(figure_path)

result = cut_polygon_by_line(polygon, line)
print(result)
plot(result)
print(result[0].intersection(result[1]))

The result is

[<shapely.geometry.polygon.Polygon object at 0x7f50dcf46d68>, 
 <shapely.geometry.polygon.Polygon object at 0x7f50dcf46da0>]
LINESTRING (-1 -1, 1 1)

enter image description here

Martin Thoma
  • 124,992
  • 159
  • 614
  • 958
  • In case you have two lines, you can just add the second line in the `linemerge()` like this: `merged = linemerge([polygon.boundary, line1, line2])` – tsveti_iko May 28 '20 at 13:01