Repository: shapely/shapely

Issue:
split ``Polygon`` by ``MultiLineString``
This code can be successfully used also for splitting `Polygon` by `MultiLineString`.

https://github.com/shapely/shapely/blob/f1b98cf1241e5d31eff6deefa1419507221d0c2f/shapely/ops.py#L365-L384

Discussion:
coveralls: ## Pull Request Test Coverage Report for [Build 12963977166](https://coveralls.io/builds/71913908)


### Details

* **2** of **3**   **(66.67%)**  changed or added relevant lines in **1** file are covered.
* No unchanged relevant lines lost coverage.
* Overall coverage increased (+**0.3%**) to **87.65%**

---

|  Changes Missing Coverage | Covered Lines | Changed/Added Lines | % |
| :-----|--------------|--------|---: |
| [shapely/ops.py](https://coveralls.io/builds/71913908/source?filename=shapely%2Fops.py#L370) | 2 | 3 | 66.67%
<!-- | **Total:** | **2** | **3** | **66.67%** | -->


|  Totals | [![Coverage Status](https://coveralls.io/builds/71913908/badge)](https://coveralls.io/builds/71913908) |
| :-- | --: |
| Change from base [Build 12798981345](https://coveralls.io/builds/71757527): |  0.3% |
| Covered Lines: | 2562 |
| Relevant Lines: | 2923 |

---
##### 💛  - [Coveralls](https://coveralls.io)


Patch:
diff --git a/CHANGES.txt b/CHANGES.txt
index b6f8dca..0a240f4 100644
--- a/CHANGES.txt
+++ b/CHANGES.txt
@@ -7,6 +7,7 @@ Changes
 Bug fixes:
 
 - Prevent crash when serializing a number > 1e100 to WKT with GEOS < 3.13. (#1907)
+- Ensure ``plot_polygon`` does not color the interiors of polygons (#1933).
 - Fixes GeoJSON serialization of empty points (#2118)
 - Fixes `__geo_interface__` handling of empty points (#2120)
 
@@ -32,7 +33,10 @@ Improvements:
 - Add option on ``invalid="fix"`` to ``from_wkb`` and ``from_wkt`` (#2094)
 - Add function ``minimum_clearance_line`` (#2106)
 - Add function ``maximum_inscribed_circle`` (#1307)
+- Add support to split polygons by multilinestrings (#2206)
 - Add function ``orient_polygons`` (#2147)
+- Add function ``coverage_simplify`` to allow topological simplification of polygonal
+  coverages (#1969)
 
 Breaking changes in GEOS 3.12:
 
diff --git a/benchmarks/benchmarks.py b/benchmarks/benchmarks.py
index 540c8ec..9c0cbd7 100644
--- a/benchmarks/benchmarks.py
+++ b/benchmarks/benchmarks.py
@@ -68,6 +68,12 @@ class ConstructorsSuite:
     def time_linearring_from_numpy(self):
         shapely.LinearRing(self.coords)
 
+    def test_polygon_from_numpy(self):
+        shapely.Polygon(self.coords)
+
+    def test_multipoint_from_numpy(self):
+        shapely.MultiPoint(self.coords)
+
 
 class ConstructiveSuite:
     """Benchmarks constructive functions on a set of 10,000 points"""
diff --git a/docs/release/1.x.rst b/docs/release/1.x.rst
index 6f673bb..356f462 100644
--- a/docs/release/1.x.rst
+++ b/docs/release/1.x.rst
@@ -361,7 +361,7 @@ Bug fixes:
 1.6.1 (2017-09-01)
 ------------------
 
-- Avoid ``STRTree`` crashes due to dangling references (#505) by maintaining
+- Avoid ``STRtree`` crashes due to dangling references (#505) by maintaining
   references to added geometries.
 - Reduce log level to debug when reporting on calls to ctypes ``CDLL()`` that
   don't succeed and are retried (#515).
diff --git a/docs/strtree.rst b/docs/strtree.rst
index 26679a8..d4803ae 100644
--- a/docs/strtree.rst
+++ b/docs/strtree.rst
@@ -1,4 +1,4 @@
-STRTree
+STRtree
 =======
 
 .. autoclass:: shapely.STRtree
diff --git a/shapely/constructive.py b/shapely/constructive.py
index 1bbbbfb..915e2ed 100644
--- a/shapely/constructive.py
+++ b/shapely/constructive.py
@@ -2,8 +2,9 @@
 
 import numpy as np
 
-from shapely import lib
+from shapely import Geometry, GeometryType, lib
 from shapely._enum import ParamEnum
+from shapely._geometry import get_parts
 from shapely.algorithms._oriented_envelope import _oriented_envelope_min_area_vectorized
 from shapely.decorators import multithreading_enabled, requires_geos
 from shapely.errors import UnsupportedGEOSVersionError
@@ -18,6 +19,7 @@ __all__ = [
     "clip_by_rect",
     "concave_hull",
     "convex_hull",
+    "coverage_simplify",
     "delaunay_triangles",
     "segmentize",
     "envelope",
@@ -1031,6 +1033,71 @@ def simplify(geometry, tolerance, preserve_topology=True, **kwargs):
         return lib.simplify(geometry, tolerance, **kwargs)
 
 
+@requires_geos("3.12.0")
+@multithreading_enabled
+def coverage_simplify(geometry, tolerance, simplify_boundary=True):
+    """Return a simplified version of an input geometry using coverage simplification.
+
+    Assumes that the geometry forms a polygonal coverage. Under this assumption, the
+    function simplifies the edges using the Visvalingam-Whyatt algorithm, while
+    preserving a valid coverage. In the most simplified case, polygons are reduced to
+    triangles.
+
+    A collection of valid polygons is considered a coverage if the polygons are:
+
+    * **Non-overlapping** - polygons do not overlap (their interiors do not intersect)
+    * **Edge-Matched** - vertices along shared edges are identical
+
+    The function allows simplification of all edges including the outer boundaries of
+    the coverage or simplification of only the inner (shared) edges.
+
+    If there are other geometry types than Polygons or MultiPolygons present, the
+    array will not undergo simplification and geometries are returned unchanged.
+
+    If the geometry is polygonal but does not form a valid coverage due to overlaps,
+    it will be simplified but it may result in invalid topology.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    tolerance : float or array_like
+        The degree of simplification roughly equal to the square root of the area
+        of triangles that will be removed.
+    simplify_boundary : bool, optional
+        By default (True), simplifies both internal edges of the coverage as well
+        as its boundary. If set to False, only simplifies internal edges.
+
+    Returns
+    -------
+    numpy.ndarray | shapely.Geometry
+
+    Examples
+    --------
+    >>> from shapely import Polygon, MultiPolygon
+    >>> poly = Polygon([(0, 0), (20, 0), (20, 10), (10, 5), (0, 10), (0, 0)])
+    >>> coverage_simplify(poly, tolerance=2)
+    <POLYGON ((0 0, 20 0, 20 10, 10 5, 0 10, 0 0))>
+    """
+    scalar = False
+    if isinstance(geometry, Geometry):
+        scalar = True
+
+    geometries = np.asarray(geometry)
+    shape = geometries.shape
+    geometries = geometries.ravel()
+
+    # create_collection acts on the inner axis
+    collections = lib.create_collection(
+        geometries, np.intc(GeometryType.GEOMETRYCOLLECTION)
+    )
+
+    simplified = lib.coverage_simplify(collections, tolerance, simplify_boundary)
+    parts = get_parts(simplified).reshape(shape)
+    if scalar:
+        return parts.item()
+    return parts
+
+
 @multithreading_enabled
 def snap(geometry, reference, tolerance, **kwargs):
     """Snap an input geometry to reference geometry's vertices.
diff --git a/shapely/geometry/multipoint.py b/shapely/geometry/multipoint.py
index 24946b2..235046a 100644
--- a/shapely/geometry/multipoint.py
+++ b/shapely/geometry/multipoint.py
@@ -1,5 +1,7 @@
 """Collections of points and related utilities."""
 
+import numpy as np
+
 import shapely
 from shapely.errors import EmptyPartError
 from shapely.geometry import point
@@ -47,16 +49,22 @@ class MultiPoint(BaseMultipartGeometry):
             return shapely.from_wkt("MULTIPOINT EMPTY")
         elif isinstance(points, MultiPoint):
             return points
+        elif len(points) == 0:
+            return shapely.from_wkt("MULTIPOINT EMPTY")
 
-        subs = []
-        for item in points:
-            p = point.Point(item)
-            if p.is_empty:
+        if isinstance(points, np.ndarray) and np.issubdtype(points.dtype, np.number):
+            subs = shapely.points(points)
+            if not subs.ndim == 1:
+                raise ValueError("Invalid values passed to MultiPoint constructor")
+            if shapely.is_empty(subs).any():
                 raise EmptyPartError("Can't create MultiPoint with empty component")
-            subs.append(p)
-
-        if len(points) == 0:
-            return shapely.from_wkt("MULTIPOINT EMPTY")
+        else:
+            subs = []
+            for item in points:
+                p = point.Point(item)
+                if p.is_empty:
+                    raise EmptyPartError("Can't create MultiPoint with empty component")
+                subs.append(p)
 
         return shapely.multipoints(subs)
 
diff --git a/shapely/geometry/polygon.py b/shapely/geometry/polygon.py
index 9347a18..489b2c0 100644
--- a/shapely/geometry/polygon.py
+++ b/shapely/geometry/polygon.py
@@ -321,6 +321,9 @@ shapely.lib.registry[3] = Polygon
 
 def orient(polygon, sign=1.0):
     """Return an oriented polygon."""
+    if polygon.is_empty:
+        return polygon
+
     s = float(sign)
     rings = []
     ring = polygon.exterior
diff --git a/shapely/ops.py b/shapely/ops.py
index 5cc0a4c..5b6a0ef 100644
--- a/shapely/ops.py
+++ b/shapely/ops.py
@@ -367,8 +367,8 @@ class SplitOp:
         """Split a Polygon with a LineString."""
         if not isinstance(poly, Polygon):
             raise GeometryTypeError("First argument must be a Polygon")
-        if not isinstance(splitter, LineString):
-            raise GeometryTypeError("Second argument must be a LineString")
+        if not isinstance(splitter, (LineString, MultiLineString)):
+            raise GeometryTypeError("Second argument must be a (Multi)LineString")
 
         union = poly.boundary.union(splitter)
 
@@ -539,7 +539,7 @@ class SplitOp:
                 )
 
         elif geom.geom_type == "Polygon":
-            if splitter.geom_type == "LineString":
+            if splitter.geom_type in ("LineString", "MultiLineString"):
                 split_func = SplitOp._split_polygon_with_line
             else:
                 raise GeometryTypeError(
diff --git a/shapely/plotting.py b/shapely/plotting.py
index 1e3070a..a9f8240 100644
--- a/shapely/plotting.py
+++ b/shapely/plotting.py
@@ -22,11 +22,14 @@ def _default_ax():
 def _path_from_polygon(polygon):
     from matplotlib.path import Path
 
+    from shapely.ops import orient
+
     if isinstance(polygon, shapely.MultiPolygon):
         return Path.make_compound_path(
             *[_path_from_polygon(poly) for poly in polygon.geoms]
         )
     else:
+        polygon = orient(polygon)
         return Path.make_compound_path(
             Path(np.asarray(polygon.exterior.coords)[:, :2]),
             *[Path(np.asarray(ring.coords)[:, :2]) for ring in polygon.interiors],
diff --git a/shapely/set_operations.py b/shapely/set_operations.py
index ebbba74..e0ed420 100644
--- a/shapely/set_operations.py
+++ b/shapely/set_operations.py
@@ -49,7 +49,7 @@ def difference(a, b, grid_size=None, **kwargs):
 
     Examples
     --------
-    >>> from shapely import box, LineString, normalize, Polygon
+    >>> from shapely import box, difference, LineString, normalize, Polygon
     >>> line = LineString([(0, 0), (2, 2)])
     >>> difference(line, LineString([(1, 1), (3, 3)]))
     <LINESTRING (0 0, 1 1)>
@@ -102,7 +102,7 @@ def intersection(a, b, grid_size=None, **kwargs):
 
     Examples
     --------
-    >>> from shapely import box, LineString, normalize, Polygon
+    >>> from shapely import box, intersection, LineString, normalize, Polygon
     >>> line = LineString([(0, 0), (2, 2)])
     >>> intersection(line, LineString([(1, 1), (3, 3)]))
     <LINESTRING (1 1, 2 2)>
@@ -150,7 +150,7 @@ def intersection_all(geometries, axis=None, **kwargs):
 
     Examples
     --------
-    >>> from shapely import LineString
+    >>> from shapely import intersection_all, LineString
     >>> line1 = LineString([(0, 0), (2, 2)])
     >>> line2 = LineString([(1, 1), (3, 3)])
     >>> intersection_all([line1, line2])
@@ -197,7 +197,7 @@ def symmetric_difference(a, b, grid_size=None, **kwargs):
 
     Examples
     --------
-    >>> from shapely import box, LineString, normalize
+    >>> from shapely import box, LineString, normalize, symmetric_difference
     >>> line = LineString([(0, 0), (2, 2)])
     >>> symmetric_difference(line, LineString([(1, 1), (3, 3)]))
     <MULTILINESTRING ((0 0, 1 1), (2 2, 3 3))>
@@ -245,7 +245,7 @@ def symmetric_difference_all(geometries, axis=None, **kwargs):
 
     Examples
     --------
-    >>> from shapely import LineString
+    >>> from shapely import LineString, symmetric_difference_all
     >>> line1 = LineString([(0, 0), (2, 2)])
     >>> line2 = LineString([(1, 1), (3, 3)])
     >>> symmetric_difference_all([line1, line2])
@@ -294,7 +294,7 @@ def union(a, b, grid_size=None, **kwargs):
 
     Examples
     --------
-    >>> from shapely import box, LineString, normalize
+    >>> from shapely import box, LineString, normalize, union
     >>> line = LineString([(0, 0), (2, 2)])
     >>> union(line, LineString([(2, 2), (3, 3)]))
     <MULTILINESTRING ((0 0, 2 2), (2 2, 3 3))>
@@ -356,7 +356,7 @@ def union_all(geometries, grid_size=None, axis=None, **kwargs):
 
     Examples
     --------
-    >>> from shapely import box, LineString, normalize, Point
+    >>> from shapely import box, LineString, normalize, Point, union_all
     >>> line1 = LineString([(0, 0), (2, 2)])
     >>> line2 = LineString([(2, 2), (3, 3)])
     >>> union_all([line1, line2])
@@ -424,7 +424,7 @@ def coverage_union(a, b, **kwargs):
 
     Examples
     --------
-    >>> from shapely import normalize, Polygon
+    >>> from shapely import coverage_union, normalize, Polygon
     >>> polygon = Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)])
     >>> normalize(coverage_union(polygon, Polygon([(1, 0), (1, 1), (2, 1), (2, 0), (1, 0)])))
     <POLYGON ((0 0, 0 1, 1 1, 2 1, 2 0, 1 0, 0 0))>
@@ -467,7 +467,7 @@ def coverage_union_all(geometries, axis=None, **kwargs):
 
     Examples
     --------
-    >>> from shapely import normalize, Polygon
+    >>> from shapely import coverage_union_all, normalize, Polygon
     >>> polygon_1 = Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)])
     >>> polygon_2 = Polygon([(1, 0), (1, 1), (2, 1), (2, 0), (1, 0)])
     >>> normalize(coverage_union_all([polygon_1, polygon_2]))
@@ -523,7 +523,7 @@ def disjoint_subset_union(a, b, **kwargs):
 
     Examples
     --------
-    >>> from shapely import normalize, Polygon
+    >>> from shapely import disjoint_subset_union, normalize, Polygon
     >>> polygon_1 = Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)])
     >>> polygon_2 = Polygon([(1, 0), (1, 1), (2, 1), (2, 0), (1, 0)])
     >>> normalize(disjoint_subset_union(polygon_1, polygon_2))
@@ -584,7 +584,7 @@ def disjoint_subset_union_all(geometries, axis=None, **kwargs):
 
     Examples
     --------
-    >>> from shapely import normalize, Polygon
+    >>> from shapely import disjoint_subset_union_all, normalize, Polygon
     >>> polygon_1 = Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)])
     >>> polygon_2 = Polygon([(1, 0), (1, 1), (2, 1), (2, 0), (1, 0)])
     >>> normalize(disjoint_subset_union_all([polygon_1, polygon_2]))
diff --git a/src/ufuncs.c b/src/ufuncs.c
index e9c7b8d..3530180 100644
--- a/src/ufuncs.c
+++ b/src/ufuncs.c
@@ -3612,6 +3612,70 @@ static PyUFuncGenericFunction to_geojson_funcs[1] = {&to_geojson_func};
 
 #endif  // GEOS_SINCE_3_10_0
 
+#if GEOS_SINCE_3_12_0
+static char coverage_simplify_dtypes[4] = {NPY_OBJECT, NPY_DOUBLE, NPY_BOOL, NPY_OBJECT};
+static void coverage_simplify_func(char** args, const npy_intp* dimensions, const npy_intp* steps,
+                                    void* data) {
+  GEOSGeometry* in1 = NULL;
+  GEOSGeometry** geom_arr;
+
+  CHECK_NO_INPLACE_OUTPUT(3);
+
+  // allocate a temporary array to store output GEOSGeometry objects
+  geom_arr = malloc(sizeof(void*) * dimensions[0]);
+  CHECK_ALLOC(geom_arr);
+
+  GEOS_INIT_THREADS;
+
+  TERNARY_LOOP {
+    CHECK_SIGNALS_THREADS(i);
+    if (errstate == PGERR_PYSIGNAL) {
+      destroy_geom_arr(ctx, geom_arr, i - 1);
+      break;
+    }
+    // get the geometry: return on error
+    if (!get_geom(*(GeometryObject**)ip1, &in1)) {
+      errstate = PGERR_NOT_A_GEOMETRY;
+      destroy_geom_arr(ctx, geom_arr, i - 1);
+      break;
+    }
+    double in2 = *(double*)ip2;
+    npy_bool in3 = !(*(npy_bool*)ip3);
+
+    int isValid = 1;
+    int numGeoms = GEOSGetNumGeometries_r(ctx, in1);
+    for (int j = 0; j < numGeoms; j++) {
+      GEOSGeometry* geom = GEOSGetGeometryN_r(ctx, in1, j);
+      if (GEOSGeomTypeId_r(ctx, geom) != GEOS_POLYGON && GEOSGeomTypeId_r(ctx, geom) != GEOS_MULTIPOLYGON) {
+          isValid = 0;
+          break;
+      }
+    }
+    if (isValid) {
+      geom_arr[i] = GEOSCoverageSimplifyVW_r(ctx, in1, in2, (int)in3);
+      if (geom_arr[i] == NULL) {
+        errstate = PGERR_GEOS_EXCEPTION;
+        destroy_geom_arr(ctx, geom_arr, i - 1);
+        break;
+      }
+    } else {
+      geom_arr[i] = GEOSGeom_clone_r(ctx, in1);
+    }
+
+  }
+
+
+  GEOS_FINISH_THREADS;
+
+  // fill the numpy array with PyObjects while holding the GIL
+  if (errstate == PGERR_SUCCESS) {
+    geom_arr_to_npy(geom_arr, args[3], steps[3], dimensions[0]);
+  }
+  free(geom_arr);
+}
+static PyUFuncGenericFunction coverage_simplify_funcs[1] = {&coverage_simplify_func};
+
+#endif  // GEOS_SINCE_3_12_0
 /*
 TODO polygonizer functions
 TODO prepared geometry predicate functions
@@ -3874,6 +3938,7 @@ int init_ufuncs(PyObject* m, PyObject* d) {
 #endif
 
 #if GEOS_SINCE_3_12_0
+  DEFINE_CUSTOM(coverage_simplify, 3);
   DEFINE_Y_Y(disjoint_subset_union);
   DEFINE_Y_b(has_m);
   DEFINE_Y_d(get_m);


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index ffd454d..67affa6 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -100,7 +100,7 @@ jobs:
         if: ${{ matrix.geos == 'main' }}
 
       - name: Set up Python ${{ matrix.python }}
-        uses: Quansight-Labs/setup-python@b9ab292c751a42bcd2bb465b7fa202ea2c3f5796  # v5.3.1
+        uses: Quansight-Labs/setup-python@869aeafb7eeb9dc48ba68acc0479e6fc3fd7ce5e  # v5.4.0
         with:
           python-version: ${{ matrix.python }}
           architecture: ${{ matrix.architecture }}
diff --git a/shapely/tests/geometry/test_multipoint.py b/shapely/tests/geometry/test_multipoint.py
index 34b43c1..aedfce5 100644
--- a/shapely/tests/geometry/test_multipoint.py
+++ b/shapely/tests/geometry/test_multipoint.py
@@ -1,6 +1,7 @@
 import numpy as np
 import pytest
 
+import shapely
 from shapely import MultiPoint, Point
 from shapely.errors import EmptyPartError
 from shapely.geometry.base import dump_coords
@@ -52,7 +53,22 @@ class TestMultiPoint(MultiGeometryTestCase):
     def test_create_multi_with_empty_component(self):
         msg = "Can't create MultiPoint with empty component"
         with pytest.raises(EmptyPartError, match=msg):
-            MultiPoint([Point(0, 0), Point()]).wkt
+            MultiPoint([Point(0, 0), Point()])
+
+        if shapely.geos_version < (3, 13, 0):
+            # for older GEOS, Point(NaN, NaN) is considered empty
+            with pytest.raises(EmptyPartError, match=msg):
+                MultiPoint([(0, 0), (np.nan, np.nan)])
+
+            with pytest.raises(EmptyPartError, match=msg):
+                MultiPoint(np.array([(0, 0), (np.nan, np.nan)]))
+
+        else:
+            result = MultiPoint([(0, 0), (np.nan, np.nan)])
+            expected = shapely.multipoints(
+                shapely.points([(0, 0), (np.nan, np.nan)], handle_nan="allow")
+            )
+            assert result == expected
 
 
 def test_multipoint_array_coercion():
diff --git a/shapely/tests/legacy/test_orient.py b/shapely/tests/legacy/test_orient.py
index 71dce26..e9df865 100644
--- a/shapely/tests/legacy/test_orient.py
+++ b/shapely/tests/legacy/test_orient.py
@@ -39,6 +39,10 @@ class OrientTestCase(unittest.TestCase):
         assert orient(linearring, 1) == linearring
         assert orient(linearring, -1) == linearring
 
+    def test_empty_polygon(self):
+        polygon = Polygon()
+        assert orient(polygon) == polygon
+
     def test_polygon(self):
         polygon = Polygon([(0, 0), (0, 1), (1, 0)])
         polygon_reversed = Polygon(polygon.exterior.coords[::-1])
diff --git a/shapely/tests/legacy/test_split.py b/shapely/tests/legacy/test_split.py
index 5ef325b..0bbf4ea 100644
--- a/shapely/tests/legacy/test_split.py
+++ b/shapely/tests/legacy/test_split.py
@@ -51,12 +51,17 @@ class TestSplitPolygon(TestSplitGeometry):
     )
 
     def test_split_poly_with_line(self):
-        # crossing at 2 points --> return 2 segments
+        # crossing at 2 points --> return 2 polygons
         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
+        # crossing twice with one linestring --> return 3 polygons
+        splitter = LineString([(1, 3), (1, -3), (1.7, -3), (1.7, 3)])
+        self.helper(self.poly_simple, splitter, 3)
+        self.helper(self.poly_hole, splitter, 3)
+
+        # 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)
@@ -71,6 +76,34 @@ class TestSplitPolygon(TestSplitGeometry):
         self.helper(self.poly_simple, splitter, 1)
         self.helper(self.poly_hole, splitter, 1)
 
+    def test_split_poly_with_multiline(self):
+        # crossing twice with a multilinestring --> return 3 polygons
+        splitter = MultiLineString([[(0.2, 3), (0.2, -3)], [(1.7, -3), (1.7, 3)]])
+        self.helper(self.poly_simple, splitter, 3)
+        self.helper(self.poly_hole, splitter, 3)
+
+        # crossing twice with a cross multilinestring --> return 4 polygons
+        splitter = MultiLineString([[(0.2, 3), (0.2, -3)], [(-3, 1), (3, 1)]])
+        self.helper(self.poly_simple, splitter, 4)
+        self.helper(self.poly_hole, splitter, 4)
+
+        # cross once, touch the boundary once --> return 2 polygons
+        splitter = MultiLineString([[(0.2, 3), (0.2, -3)], [(0, 2), (5, 2)]])
+        self.helper(self.poly_simple, splitter, 2)
+        self.helper(self.poly_hole, splitter, 2)
+
+        # cross once, inside the polygon once --> return 2 polygons
+        splitter = MultiLineString(
+            [[(0.2, 3), (0.2, -3)], [(1.2, 1.2), (1.7, 1.7), (3, 2)]]
+        )
+        self.helper(self.poly_simple, splitter, 2)
+        self.helper(self.poly_hole, splitter, 2)
+
+        # cross once, outside the polygon once --> return 2 polygons
+        splitter = MultiLineString([[(0.2, 3), (0.2, -3)], [(0, 3), (3, 3), (3, 0)]])
+        self.helper(self.poly_simple, splitter, 2)
+        self.helper(self.poly_hole, splitter, 2)
+
     def test_split_poly_with_other(self):
         with pytest.raises(GeometryTypeError):
             split(self.poly_simple, Point(1, 1))
diff --git a/shapely/tests/test_constructive.py b/shapely/tests/test_constructive.py
index 352f6ef..0c6b5e2 100644
--- a/shapely/tests/test_constructive.py
+++ b/shapely/tests/test_constructive.py
@@ -1097,6 +1097,86 @@ def test_concave_hull_kwargs():
     assert shapely.get_num_coordinates(result4) < shapely.get_num_coordinates(result3)
 
 
+@pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12")
+@pytest.mark.parametrize("geometry", all_types)
+def test_coverage_simplify_scalars(geometry):
+    actual = shapely.coverage_simplify(geometry, 0.0)
+    assert isinstance(actual, Geometry)
+    assert shapely.get_type_id(actual) == shapely.get_type_id(geometry)
+    # Anything other than MultiPolygon or a GeometryCollection is returned as-is
+    if shapely.get_type_id(geometry) not in (3, 6):
+        assert actual.equals(geometry)
+
+
+@pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12")
+@pytest.mark.parametrize("geometry", all_types)
+def test_coverage_simplify_geom_types(geometry):
+    actual = shapely.coverage_simplify([geometry, geometry], 0.0)
+    assert isinstance(actual, np.ndarray)
+    assert actual.shape == (2,)
+    assert (shapely.get_type_id(actual) == shapely.get_type_id(geometry)).all()
+
+
+@pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12")
+def test_coverage_simplify_multipolygon():
+    mp = MultiPolygon(
+        [
+            Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)]),
+            Polygon([(2, 2), (2, 3), (3, 3), (3, 2), (2, 2)]),
+        ]
+    )
+    actual = shapely.coverage_simplify(mp, 1)
+    assert actual.equals(
+        shapely.from_wkt(
+            "MULTIPOLYGON (((0 1, 1 1, 1 0, 0 1)), ((2 3, 3 3, 3 2, 2 3)))"
+        )
+    )
+
+
+@pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12")
+def test_coverage_simplify_array():
+    polygons = np.array(
+        [
+            shapely.Polygon([(0, 0), (20, 0), (20, 10), (10, 5), (0, 10), (0, 0)]),
+            shapely.Polygon([(0, 10), (10, 5), (20, 10), (20, 20), (0, 20), (0, 10)]),
+        ]
+    )
+    low_tolerance = shapely.coverage_simplify(polygons, 1)
+    mid_tolerance = shapely.coverage_simplify(polygons, 8)
+    high_tolerance = shapely.coverage_simplify(polygons, 10)
+
+    assert shapely.equals(low_tolerance, shapely.normalize(polygons)).all()
+    assert shapely.equals(
+        mid_tolerance,
+        shapely.from_wkt(
+            [
+                "POLYGON ((20 10, 0 10, 0 0, 20 0, 20 10))",
+                "POLYGON ((20 10, 0 10, 0 20, 20 20, 20 10))",
+            ]
+        ),
+    ).all()
+    assert shapely.equals(
+        high_tolerance,
+        shapely.from_wkt(
+            [
+                "POLYGON ((20 10, 0 10, 20 0, 20 10))",
+                "POLYGON ((20 10, 0 10, 0 20, 20 10))",
+            ]
+        ),
+    ).all()
+
+    no_boundary = shapely.coverage_simplify(polygons, 10, simplify_boundary=False)
+    assert shapely.equals(
+        no_boundary,
+        shapely.from_wkt(
+            [
+                "POLYGON ((20 10, 0 10, 0 0, 20 0, 20 10))",
+                "POLYGON ((20 10, 0 10, 0 20, 20 20, 20 10))",
+            ]
+        ),
+    ).all()
+
+
 @pytest.mark.skipif(shapely.geos_version < (3, 12, 0), reason="GEOS < 3.12")
 def test_voronoi_polygons_ordered():
     mp = MultiPoint([(3.0, 1.0), (3.0, 2.0), (1.0, 2.0), (1.0, 1.0)])
diff --git a/shapely/tests/test_plotting.py b/shapely/tests/test_plotting.py
index b2b9c2a..78a4175 100644
--- a/shapely/tests/test_plotting.py
+++ b/shapely/tests/test_plotting.py
@@ -1,7 +1,15 @@
 import pytest
 from numpy.testing import assert_allclose
 
-from shapely import LineString, MultiLineString, Point, box, get_coordinates
+from shapely import (
+    LineString,
+    MultiLineString,
+    MultiPolygon,
+    Point,
+    box,
+    get_coordinates,
+)
+from shapely.ops import orient
 from shapely.plotting import patch_from_polygon, plot_line, plot_points, plot_polygon
 
 pytest.importorskip("matplotlib")
@@ -48,7 +56,7 @@ def test_plot_polygon_with_interior():
     poly = box(0, 0, 1, 1).difference(box(0.2, 0.2, 0.5, 0.5))
     artist, _ = plot_polygon(poly)
     plot_coords = artist.get_path().vertices
-    assert_allclose(plot_coords, get_coordinates(poly))
+    assert_allclose(plot_coords, get_coordinates(orient(poly)))
 
 
 def test_plot_multipolygon():
@@ -58,6 +66,16 @@ def test_plot_multipolygon():
     assert_allclose(plot_coords, get_coordinates(poly))
 
 
+def test_plot_multipolygon_with_interior():
+    poly1 = box(0, 0, 1, 1).difference(box(0.2, 0.2, 0.5, 0.5))
+    poly2 = box(3, 3, 6, 6).difference(box(4, 4, 5, 5))
+    poly = MultiPolygon([poly1, poly2])
+
+    artist, _ = plot_polygon(poly)
+    plot_coords = artist.get_path().vertices
+    assert_allclose(plot_coords, get_coordinates(orient(poly)))
+
+
 def test_plot_line():
     line = LineString([(0, 0), (1, 0), (1, 1)])
     artist, _ = plot_line(line)
diff --git a/shapely/tests/test_strtree.py b/shapely/tests/test_strtree.py
index d54d4c5..9911917 100644
--- a/shapely/tests/test_strtree.py
+++ b/shapely/tests/test_strtree.py
@@ -118,9 +118,6 @@ def test_geometries_property():
     assert_geometries_equal(point, tree.geometries[0])
 
 
-# TODO(shapely-2.0) this fails on Appveyor, see
-# https://github.com/shapely/shapely/pull/983#issuecomment-718557666
-@pytest.mark.skipif(sys.platform.startswith("win32"), reason="does not run on Appveyor")
 def test_pickle_persistence(tmp_path):
     # write the pickeled tree to another process; the process should not crash
     tree = STRtree([Point(i, i).buffer(0.1) for i in range(3)])