Repository: Toblerity/Shapely
Pull request #371 (2016-03-09 12:32:49)

Add shapely.ops.split function.
This branch provides the functionality for splitting a geometry by another geometry as requested on #293,
based on the implementation of ST_Split used in PostGIS. 
I added the function **shapely.ops.split**, some tests and the documentation.


Patch:
diff --git a/docs/code/minimum_rotated_rectangle.py b/docs/code/minimum_rotated_rectangle.py
deleted file mode 100644
index 7eb672c..0000000
--- a/docs/code/minimum_rotated_rectangle.py
+++ /dev/null
@@ -1,49 +0,0 @@
-from shapely.geometry import MultiPoint, Polygon, LineString
-import matplotlib.pyplot as plt
-from descartes.patch import PolygonPatch
-
-from figures import SIZE
-
-fig = plt.figure(1, figsize=SIZE, dpi=90)
-fig.set_frameon(True)
-
-# 1
-ax = fig.add_subplot(121)
-
-mp = MultiPoint([(0, 0), (0.5, 1.5), (1, 0.5), (0.5, 0.5)])
-rect = mp.minimum_rotated_rectangle
-
-for p in mp:
-	ax.plot(p.x, p.y, 'o', color='#999999')
-patch = PolygonPatch(rect, facecolor='#6699cc', edgecolor='#6699cc', alpha=0.5, zorder=2)
-ax.add_patch(patch)
-ax.set_title('a) MultiPoint')
-
-xr = [-1, 2]
-yr = [-1, 2]
-ax.set_xlim(*xr)
-ax.set_xticks(range(*xr) + [xr[-1]])
-ax.set_ylim(*yr)
-ax.set_yticks(range(*yr) + [yr[-1]])
-ax.set_aspect(1)
-
-# 2
-ax = fig.add_subplot(122)
-ls = LineString([(-0.5, 1.2), (0.5, 0), (1, 1), (1.5, 0), (1.5, 0.5)])
-rect = ls.minimum_rotated_rectangle
-
-ax.plot(*ls.xy, color='#333333', linewidth=3, alpha=0.5, zorder=2)
-patch = PolygonPatch(rect, facecolor='#6699cc', edgecolor='#6699cc', alpha=0.5, zorder=2)
-ax.add_patch(patch)
-
-xr = [-1, 2]
-yr = [-1, 2]
-ax.set_xlim(*xr)
-ax.set_xticks(range(*xr) + [xr[-1]])
-ax.set_ylim(*yr)
-ax.set_yticks(range(*yr) + [yr[-1]])
-ax.set_aspect(1)
-
-ax.set_title('b) LineString')
-
-plt.show()
diff --git a/docs/manual.rst b/docs/manual.rst
index e540011..0508f8e 100644
--- a/docs/manual.rst
+++ b/docs/manual.rst
@@ -1576,27 +1576,6 @@ Figure 10. Convex hull (blue) of 2 points (left) and of 6 points (right).
   >>> MultiPoint([(0, 0), (1, 1)]).envelope
   <shapely.geometry.polygon.Polygon object at 0x...>
 
-.. attribute:: object.minimum_rotated_rectangle
-  
-  Returns the general minimum bounding rectangle that contains the object. 
-  Unlike envelope this rectangle is not constrained to be parallel to the 
-  coordinate axes. If the convex hull of the object is a degenerate (line or point) 
-  this degenerate is returned.
-
-  `New in Shapely 1.6.0`
-
-.. code-block:: pycon
-  
-  >>> Point(0, 0).minimum_rotated_rectangle
-  <shapely.geometry.point.Point object at 0x...>
-  >>> MultiPoint([(0,0),(1,1),(2,0.5)]).minimum_rotated_rectangle
-  <shapely.geometry.polygon.Polygon object at 0x...>
-
-.. plot:: code/minimum_rotated_rectangle.py
-
-Figure 11. Minimum rotated rectangle for a multipoint feature (left) and a 
-linestring feature (right).
-
 .. method:: object.parallel_offset(distance, side, resolution=16, join_style=1, mitre_limit=5.0)
 
   Returns a LineString or MultiLineString geometry at a distance from the
@@ -1621,7 +1600,7 @@ linestring feature (right).
 
 .. plot:: code/parallel_offset.py
 
-Figure 12. Three styles of parallel offset lines on the left side of a simple
+Figure 11. Three styles of parallel offset lines on the left side of a simple
 line string (its starting point shown as a circle) and one offset on the right
 side, a multipart.
 
@@ -1629,7 +1608,7 @@ The effect of the `mitre_limit` parameter is shown below.
 
 .. plot:: code/parallel_offset_mitre.py
 
-Figure 13. Large and small mitre_limit values for left and right offsets.
+Figure 12. Large and small mitre_limit values for left and right offsets.
 
 .. method:: object.simplify(tolerance, preserve_topology=True)
 
@@ -1656,7 +1635,7 @@ Douglas-Peucker algorithm [6]_ is used.
 
 .. plot:: code/simplify.py
 
-Figure 14. Simplification of a nearly circular polygon using a tolerance of 0.2
+Figure 13. Simplification of a nearly circular polygon using a tolerance of 0.2
 (left) and 0.5 (right).
 
 .. note::
@@ -1764,7 +1743,7 @@ preserved or supported by 3D affine transformations.
 
   .. plot:: code/rotate.py
 
-  Figure 15. Rotation of a `LineString` (gray) by an angle of 90°
+  Figure 14. Rotation of a `LineString` (gray) by an angle of 90°
   counter-clockwise (blue) using different origins.
 
 .. function:: shapely.affinity.scale(geom, xfact=1.0, yfact=1.0, zfact=1.0, origin='center')
@@ -1806,7 +1785,7 @@ preserved or supported by 3D affine transformations.
 
   .. plot:: code/scale.py
 
-  Figure 16. Scaling of a gray triangle to blue result: a) by a factor of 1.5
+  Figure 15. Scaling of a gray triangle to blue result: a) by a factor of 1.5
   along x-direction, with reflection across y-axis; b) by a factor of 2 along
   x-direction with custom origin at (1, 1).
 
@@ -1838,7 +1817,7 @@ preserved or supported by 3D affine transformations.
 
   .. plot:: code/skew.py
 
-  Figure 17. Skewing of a gray "R" to blue result: a) by a shear angle of 20°
+  Figure 16. Skewing of a gray "R" to blue result: a) by a shear angle of 20°
   along the x-direction and an origin at (1, 1); b) by a shear angle of 30°
   along the y-direction, using default origin.
 
@@ -2147,6 +2126,36 @@ paths between two lineal geometries.
   >>> backward.wkt
   'MULTILINESTRING ((10 5, 20 5))'
 
+Splitting
+---------
+
+The :func:`~shapely.ops.split` function in `shapely.ops` splits a geometry by another geometry.
+
+.. function:: shapely.ops.split(geom, splitter)
+
+   Splits a geometry by another geometry and returns a collection of geometries. This function is the theoretical
+   opposite of the union of the split geometry parts. If the splitter does not split the geometry, a  collection with a single geometry equal to the input geometry is returned.
+  
+   The function supports:
+
+   * Splitting a (Multi)LineString by a (Multi)Point or (Multi)LineString or (Multi)Polygon boundary
+   
+   * Splitting a (Multi)Polygon by a LineString
+
+   It may be convenient to snap the splitter with low tolerance to the geometry. For example in the case of splitting a line by a point, the point must be exactly on the line, for the line to be correctly split.
+   When splitting a line by a polygon, the boundary of the polygon is used for the operation.
+   When splitting a line by another line, a ValueError is raised if the two overlap at some segment.
+
+   `New in version 1.6.0`
+
+.. code-block:: pycon
+
+  >>> pt = Point((1, 1))
+  >>> line = LineString([(0,0), (2,2)])
+  >>> result = split(line, pt)
+  >>> result.wkt
+  'GEOMETRYCOLLECTION (LINESTRING (0 0, 1 1), LINESTRING (1 1, 2 2))'
+
 Prepared Geometry Operations
 ----------------------------
 
diff --git a/shapely/ops.py b/shapely/ops.py
index c27b9d4..2c7a79b 100644
--- a/shapely/ops.py
+++ b/shapely/ops.py
@@ -12,11 +12,11 @@ from ctypes import byref, c_void_p, c_double
 
 from shapely.geos import lgeos
 from shapely.geometry.base import geom_factory, BaseGeometry
-from shapely.geometry import asShape, asLineString, asMultiLineString, Point, \
-                             LineString
+from shapely.geometry import asShape, asLineString, asMultiLineString, Point, MultiPoint,\
+                             LineString, MultiLineString, Polygon, MultiPolygon, GeometryCollection
 
 __all__ = ['cascaded_union', 'linemerge', 'operator', 'polygonize',
-           'polygonize_full', 'transform', 'unary_union', 'triangulate']
+           'polygonize_full', 'transform', 'unary_union', 'triangulate', 'split']
 
 
 class CollectionOperator(object):
@@ -324,3 +324,151 @@ def shared_paths(g1, g2):
     if not isinstance(g2, LineString):
         raise TypeError("Second geometry must be a LineString")
     return(geom_factory(lgeos.methods['shared_paths'](g1._geom, g2._geom)))
+
+
+class SplitOp(object):
+
+    @staticmethod
+    def _split_polygon_with_line(poly, splitter):
+        """Split a Polygon with a LineString"""
+
+        assert(isinstance(poly, Polygon))
+        assert(isinstance(splitter, LineString))
+
+        union = poly.boundary.union(splitter)
+
+        # some polygonized geometries may be holes, we do not want them
+        # that's why we test if the original polygon (poly) contains 
+        # an inner point of polygonized geometry (pg)
+        return [pg for pg in polygonize(union) if poly.contains(pg.representative_point())]
+
+    @staticmethod
+    def _split_line_with_line(line, splitter):
+        """Split a LineString with another (Multi)LineString or (Multi)Polygon"""
+        
+        # if splitter is a polygon, pick it's boundary
+        if splitter.type in ('Polygon', 'MultiPolygon'):
+            splitter = splitter.boundary
+
+        assert(isinstance(line, LineString))
+        assert(isinstance(splitter, LineString) or isinstance(splitter, MultiLineString))
+        
+        if splitter.crosses(line):
+            # The lines cross --> return multilinestring from the split
+            return line.difference(splitter)
+        elif splitter.relate_pattern(line, '1********'):
+            # The lines overlap at some segment (linear intersection of interiors)
+            raise ValueError('Input geometry segment overlaps with the splitter.')
+        else:
+            # The lines do not cross --> return collection with identity line
+            return [line]
+
+    @staticmethod  
+    def _split_line_with_point(line, splitter):
+        """Split a LineString with a Point"""
+
+        assert(isinstance(line, LineString))
+        assert(isinstance(splitter, Point))
+
+        # check if point is in the interior of the line
+        if not line.relate_pattern(splitter, '0********'):
+            # point not on line interior --> return collection with single identity line
+            # (REASONING: Returning a list with the input line reference and creating a 
+            # GeometryCollection at the general split function prevents unnecessary copying 
+            # of linestrings in multipoint splitting function)
+            return [line]
+
+        # point is on line, get the distance from the first point on line
+        distance_on_line = line.project(splitter)
+        coords = list(line.coords)
+        # split the line at the point and create two new lines
+        # TODO: can optimize this by accumulating the computed point-to-point distances
+        for i, p in enumerate(coords):
+            pd = line.project(Point(p))
+            if pd == distance_on_line:
+                return [
+                    LineString(coords[:i+1]), 
+                    LineString(coords[i:])
+                ]
+            elif distance_on_line < pd:
+                # we must interpolate here because the line might use 3D points
+                cp = line.interpolate(distance_on_line)
+                ls1_coords = coords[:i]
+                ls1_coords.append(cp.coords[0])
+                ls2_coords = [cp.coords[0]]
+                ls2_coords.extend(coords[i:])
+                return [LineString(ls1_coords), LineString(ls2_coords)]
+
+    @staticmethod
+    def _split_line_with_multipoint(line, splitter):
+        """Split a LineString with a MultiPoint"""
+
+        assert(isinstance(line, LineString))
+        assert(isinstance(splitter, MultiPoint))
+        
+        chunks = [line]
+        for pt in splitter.geoms:
+            new_chunks = []
+            for chunk in chunks:
+                # add the newly split 2 lines or the same line if not split
+                new_chunks.extend(SplitOp._split_line_with_point(chunk, pt))
+            chunks = new_chunks
+        
+        return chunks
+    
+    @staticmethod
+    def split(geom, splitter):
+        """
+        Splits a geometry by another geometry and returns a collection of geometries. This function is the theoretical
+        opposite of the union of the split geometry parts. If the splitter does not split the geometry, a collection
+        with a single geometry equal to the input geometry is returned.
+        The function supports:
+          - Splitting a (Multi)LineString by a (Multi)Point or (Multi)LineString or (Multi)Polygon
+          - Splitting a (Multi)Polygon by a LineString
+
+        It may be convenient to snap the splitter with low tolerance to the geometry. For example in the case
+        of splitting a line by a point, the point must be exactly on the line, for the line to be correctly split.
+        When splitting a line by a polygon, the boundary of the polygon is used for the operation.
+        When splitting a line by another line, a ValueError is raised if the two overlap at some segment.
+
+        Parameters
+        ----------
+        geom : geometry
+            The geometry to be split
+        splitter : geometry
+            The geometry that will split the input geom
+
+        Example
+        -------
+        >>> pt = Point((1, 1))
+        >>> line = LineString([(0,0), (2,2)])
+        >>> result = split(line, pt)
+        >>> result.wkt
+        'GEOMETRYCOLLECTION (LINESTRING (0 0, 1 1), LINESTRING (1 1, 2 2))'
+        """
+
+        if geom.type in ('MultiLineString', 'MultiPolygon'):
+             return GeometryCollection([i for part in geom.geoms for i in SplitOp.split(part, splitter).geoms])
+
+        elif geom.type == 'LineString':
+            if splitter.type in ('LineString', 'MultiLineString', 'Polygon', 'MultiPolygon'):
+                split_func = SplitOp._split_line_with_line
+            elif splitter.type in ('Point'):
+                split_func = SplitOp._split_line_with_point
+            elif splitter.type in ('MultiPoint'):
+                split_func =  SplitOp._split_line_with_multipoint
+            else:
+                raise ValueError("Splitting a LineString with a %s is not supported" % splitter.type)
+
+        elif geom.type == 'Polygon':
+            if splitter.type == 'LineString':
+                split_func = SplitOp._split_polygon_with_line
+            else:
+                raise ValueError("Splitting a Polygon with a %s is not supported" % splitter.type)
+
+        else:
+            raise ValueError("Splitting %s geometry is not supported" % geom.type)
+
+        return GeometryCollection(split_func(geom, splitter))
+
+split = SplitOp.split
\ No newline at end of file


Test patch:
diff --git a/tests/test_split.py b/tests/test_split.py
new file mode 100644
index 0000000..dcaae0b
--- /dev/null
+++ b/tests/test_split.py
@@ -0,0 +1,192 @@
+from shapely.ops import split
+
+from . import unittest
+
+from shapely.geometry import Point, MultiPoint, LineString, MultiLineString, Polygon, MultiPolygon, GeometryCollection
+from shapely.ops import cascaded_union, linemerge
+
+class TestSplitGeometry(unittest.TestCase):
+	# helper class for testing below
+	def helper(self, geom, splitter, expected_chunks):	
+		s = split(geom, splitter)
+		self.assertEqual(s.type, "GeometryCollection")
+		self.assertEqual(len(s), expected_chunks)
+		if expected_chunks > 1:
+			# split --> expected collection that when merged is again equal to original geometry
+			if s.geoms[0].type == 'LineString':
+				self.assertTrue(linemerge(s).simplify(0.000001).equals(geom))
+			elif s.geoms[0].type == 'Polygon':
+				union = cascaded_union(s).simplify(0.000001)
+				self.assertTrue(union.equals(geom))
+				self.assertEqual(union.area, geom.area)
+			else:
+				raise ValueError
+		elif expected_chunks == 1:
+			# not split --> expected equal to line
+			self.assertTrue(s[0].equals(geom))
+
+class TestSplitPolygon(TestSplitGeometry):
+	poly_simple = Polygon([(0, 0), (2, 0), (2, 2), (0, 2), (0, 0)])
+	poly_hole = Polygon([(0, 0), (2, 0), (2, 2), (0, 2), (0, 0)], [[(0.5, 0.5), (0.5, 1.5), (1.5, 1.5), (1.5, 0.5), (0.5, 0.5)]])
+
+	def test_split_poly_with_line(self):
+		# crossing at 2 points --> return 2 segments
+		splitter = LineString([(1, 3), (1, -3)])
+		self.helper(self.poly_simple, splitter, 2)
+		self.helper(self.poly_hole, splitter, 2)
+
+		# touching the boundary--> return equal
+		splitter = LineString([(0, 2), (5, 2)])
+		self.helper(self.poly_simple, splitter, 1)
+		self.helper(self.poly_hole, splitter, 1)
+
+		# inside the polygon --> return equal
+		splitter = LineString([(0.2, 0.2), (1.7, 1.7), (3, 2)])
+		self.helper(self.poly_simple, splitter, 1)
+		self.helper(self.poly_hole, splitter, 1)
+
+		# outside the polygon --> return equal
+		splitter = LineString([(0, 3), (3, 3) , (3, 0)])
+		self.helper(self.poly_simple, splitter, 1)
+		self.helper(self.poly_hole, splitter, 1)
+
+	def test_split_poly_with_other(self):
+		with self.assertRaises(ValueError):
+			split(self.poly_simple, Point(1, 1))
+		with self.assertRaises(ValueError):
+			split(self.poly_simple, MultiPoint([(1, 1), (3, 4)]))
+		with self.assertRaises(ValueError):
+			split(self.poly_simple, self.poly_hole)
+
+class TestSplitLine(TestSplitGeometry):
+	ls = LineString([(0, 0), (1.5, 1.5), (3.0, 4.0)])
+
+	def test_split_line_with_point(self):
+		# point on line interior --> return 2 segments
+		splitter = Point(1, 1)
+		self.helper(self.ls, splitter, 2)
+
+		# point on line point --> return 2 segments
+		splitter = Point(1.5, 1.5)
+		self.helper(self.ls, splitter, 2)
+
+		# point on boundary --> return equal
+		splitter = Point(3, 4)
+		self.helper(self.ls, splitter, 1)
+
+		# point on exterior of line --> return equal
+		splitter = Point(2, 2)
+		self.helper(self.ls, splitter, 1)
+
+	def test_split_line_with_multipoint(self):
+		# points on line interior --> return 4 segments
+		splitter = MultiPoint([(1,1), (1.5, 1.5), (0.5, 0.5)])
+		self.helper(self.ls, splitter, 4)
+
+		# points on line interior and boundary -> return 2 segments
+		splitter = MultiPoint([(1, 1), (3, 4)])
+		self.helper(self.ls, splitter, 2)
+
+		# point on linear interior but twice --> return 2 segments
+		splitter = MultiPoint([(1, 1), (1.5, 1.5), (1, 1)])
+		self.helper(self.ls, splitter, 3)
+
+	def test_split_line_with_line(self):
+		# crosses at one point --> return 2 segments
+		splitter = LineString([(0, 1), (1, 0)])
+		self.helper(self.ls, splitter, 2)	
+
+		# crosses at two points --> return 3 segments
+		splitter = LineString([(0, 1), (1, 0), (1, 2)])
+		self.helper(self.ls, splitter, 3)
+
+		# overlaps --> raise
+		splitter = LineString([(0, 0), (15, 15)])
+		with self.assertRaises(ValueError):
+			self.helper(self.ls, splitter, 1)
+
+		# does not cross --> return equal
+		splitter = LineString([(0, 1), (0, 2)])
+		self.helper(self.ls, splitter, 1)
+
+		# is touching the boundary --> return equal
+		splitter = LineString([(-1, 1), (1, -1)])
+		self.assertTrue(splitter.touches(self.ls))
+		self.helper(self.ls, splitter, 1)
+
+	def test_split_line_with_multiline(self):
+		# crosses at one point --> return 2 segments
+		splitter = MultiLineString([[(0, 1), (1, 0)], [(0, 0), (2, -2)]])
+		self.helper(self.ls, splitter, 2)	
+
+		# crosses at two points --> return 3 segments
+		splitter = MultiLineString([[(0, 1), (1, 0)], [(0, 2), (2, 0)]])
+		self.helper(self.ls, splitter, 3)
+
+		# crosses at three points --> return 4 segments
+		splitter = MultiLineString([[(0, 1), (1, 0)], [(0, 2), (2, 0), (2.2, 3.2)]])
+		self.helper(self.ls, splitter, 4)
+
+		# overlaps --> raise
+		splitter = MultiLineString([[(0, 0), (1.5, 1.5)], [(1.5, 1.5), (3, 4)]])
+		with self.assertRaises(ValueError):
+			self.helper(self.ls, splitter, 1)
+
+		# does not cross --> return equal
+		splitter = MultiLineString([[(0, 1), (0, 2)], [(1, 0), (2, 0)]])
+		self.helper(self.ls, splitter, 1)
+
+	def test_split_line_with_polygon(self):
+		# crosses at two points --> return 3 segments
+		splitter = Polygon([(1, 0), (1, 2), (2, 2), (2, 0), (1, 0)])
+		self.helper(self.ls, splitter, 3)
+
+		# crosses at one point and touches boundary --> return 2 segments
+		splitter = Polygon([(0, 0), (1, 2), (2, 2), (1, 0), (0, 0)])
+		self.helper(self.ls, splitter, 2)
+
+		# exterior crosses at one point and touches at (0, 0)
+		# interior crosses at two points
+		splitter = Polygon([(0, 0), (2, 0), (2, 2), (0, 2), (0, 0)], [[(0.5, 0.5), (0.5, 1.5), (1.5, 1.5), (1.5, 0.5), (0.5, 0.5)]])
+		self.helper(self.ls, splitter, 4)
+
+	def test_split_line_with_multipolygon(self):
+		poly1 = Polygon([(0, 0), (2, 0), (2, 2), (0, 2), (0, 0)]) # crosses at one point and touches at (0, 0)
+		poly2 = Polygon([(0.5, 0.5), (0.5, 1.5), (1.5, 1.5), (1.5, 0.5), (0.5, 0.5)]) # crosses at two points
+		poly3 = Polygon([(0, 0), (0, -2), (-2, -2), (-2, 0), (0, 0)]) # not crossing
+		splitter = MultiPolygon([poly1, poly2, poly3])
+		self.helper(self.ls, splitter, 4)
+
+class TestSplitMulti(TestSplitGeometry):
+	
+	def test_split_multiline_with_point(self):
+		# a cross-like multilinestring with a point in the middle --> return 4 line segments
+		l1 = LineString([(0, 1), (2, 1)])
+		l2 = LineString([(1, 0), (1, 2)])
+		ml = MultiLineString([l1, l2])
+		splitter = Point((1, 1))
+		self.helper(ml, splitter, 4)
+		
+	def test_split_multiline_with_multipoint(self):
+		# a cross-like multilinestring with a point in middle, a point on one of the lines and a point in the exterior 
+		# --> return 4+1 line segments
+		l1 = LineString([(0, 1), (3, 1)])
+		l2 = LineString([(1, 0), (1, 2)])
+		ml = MultiLineString([l1, l2])
+		splitter = MultiPoint([(1, 1), (2, 1), (4, 2)])
+		self.helper(ml, splitter, 5)
+				
+	def test_split_multipolygon_with_line(self):
+		# two polygons with a crossing line --> return 4 triangles
+		poly1 = Polygon([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]) 
+		poly2 = Polygon([(1, 1), (1, 2), (2, 2), (2, 1), (1, 1)])
+		mpoly = MultiPolygon([poly1, poly2])
+		ls = LineString([(-1, -1), (3, 3)])
+		self.helper(mpoly, ls, 4)
+
+		# two polygons away from the crossing line --> return identity
+		poly1 = Polygon([(10, 10), (10, 11), (11, 11), (11, 10), (10, 10)]) 
+		poly2 = Polygon([(-10, -10), (-10, -11), (-11, -11), (-11, -10), (-10, -10)])
+		mpoly = MultiPolygon([poly1, poly2])
+		ls = LineString([(-1, -1), (3, 3)])
+		self.helper(mpoly, ls, 2)
\ No newline at end of file