FIX: create a single inverted polygon when no exteriors found

Fixes #1149, #1674, #2160,

Author: rcomer

lib/cartopy/crs.py

@@ -21,6 +21,7 @@
 import pyproj
 from pyproj import Transformer
 from pyproj.exceptions import ProjError
+import shapely
 import shapely.geometry as sgeom
 from shapely.prepared import prep
 
@@ -1214,34 +1215,52 @@ def _rings_to_multi_polygon(self, rings, is_ccw):
             y3 -= by
             x4 += bx
             y4 += by
+
+            interior_polys = []
+
             for ring in interior_rings:
-                # Use shapely buffer in an attempt to fix invalid geometries
-                polygon = sgeom.Polygon(ring).buffer(0)
-                if not polygon.is_empty and polygon.is_valid:
+                polygon = shapely.make_valid(sgeom.Polygon(ring))
+                if not polygon.is_empty:
+                    if isinstance(polygon, sgeom.Polygon):
+                        interior_polys.append(polygon)
+                    elif isinstance(polygon, sgeom.MultiPolygon):
+                        interior_polys.extend(polygon.geoms)
+                    elif isinstance(polygon, sgeom.GeometryCollection):
+                        for geom in polygon.geoms:
+                            if isinstance(geom, sgeom.Polygon):
+                                interior_polys.append(geom)
+                    else:
+                        # make_valid may produce some linestrings.  Ignore these
+                        continue
+
                     x1, y1, x2, y2 = polygon.bounds
                     bx = (x2 - x1) * 0.1
                     by = (y2 - y1) * 0.1
                     x1 -= bx
                     y1 -= by
                     x2 += bx
                     y2 += by
-                    box = sgeom.box(min(x1, x3), min(y1, y3),
-                                    max(x2, x4), max(y2, y4))
 
-                    # Invert the polygon
-                    polygon = box.difference(polygon)
+                    x3 = min(x1, x3)
+                    x4 = max(x2, x4)
+                    y3 = min(y1, y3)
+                    y4 = max(y2, y4)
 
-                    # Intersect the inverted polygon with the boundary
-                    polygon = boundary_poly.intersection(polygon)
+            box = sgeom.box(x3, y3, x4, y4, ccw=is_ccw)
 
-                    if not polygon.is_empty:
-                        polygon_bits.append(polygon)
+            # Invert the polygons
+            polygon = box.difference(sgeom.MultiPolygon(interior_polys))
 
-        if polygon_bits:
-            multi_poly = sgeom.MultiPolygon(polygon_bits)
-        else:
-            multi_poly = sgeom.MultiPolygon()
-        return multi_poly
+            # Intersect the inverted polygon with the boundary
+            polygon = boundary_poly.intersection(polygon)
+
+            if not polygon.is_empty:
+                if isinstance(polygon, sgeom.MultiPolygon):
+                    polygon_bits.extend(polygon.geoms)
+                else:
+                    polygon_bits.append(polygon)
+
+        return sgeom.MultiPolygon(polygon_bits)
 
     def quick_vertices_transform(self, vertices, src_crs):
         """

lib/cartopy/tests/mpl/baseline_images/mpl/test_features/ocean_lambertazimuthalequalarea.png

@@ -46,6 +46,18 @@ def test_natural_earth_custom():
     return ax.figure
 
 
+@pytest.mark.filterwarnings("ignore:Downloading")
+@pytest.mark.natural_earth
+@pytest.mark.mpl_image_compare(filename='ocean_lambertazimuthalequalarea.png')
+def test_ocean_lambertazimuthalequalarea():
+    proj = ccrs.LambertAzimuthalEqualArea(central_latitude=45, central_longitude=-100)
+    ax = plt.axes(projection=proj)
+    ax.add_feature(cfeature.OCEAN)
+    ax.coastlines()
+    ax.set_extent([-140, -70, 20, 60])
+    return ax.figure
+
+
 @pytest.mark.network
 @pytest.mark.skipif(not _HAS_PYKDTREE_OR_SCIPY, reason='pykdtree or scipy is required')
 @pytest.mark.mpl_image_compare(filename='gshhs_coastlines.png', tolerance=0.95)

lib/cartopy/tests/mpl/test_features.py

@@ -446,6 +446,18 @@ def test_inverted_poly_simple_hole(self):
         self._assert_bounds(polygon.interiors[0].bounds,
                             - 1.2e7, -1.2e7, 1.2e7, 1.2e7, 1e6)
 
+    def test_inverted_poly_multi_hole(self):
+        # Adapted from 1149
+        proj = ccrs.LambertAzimuthalEqualArea(
+            central_latitude=45, central_longitude=-100)
+        poly = sgeom.Polygon([(-180, -80), (180, -80), (180, 90), (-180, 90)],
+                             [[(-50, -50), (-50, 0), (0, 0), (0, -50)]])
+        multi_polygon = proj.project_geometry(poly)
+        # Should project to single polygon with multiple holes
+        assert len(multi_polygon.geoms) == 1
+        assert len(multi_polygon.geoms[0].interiors) >= 2
+
+
     def test_inverted_poly_clipped_hole(self):
         proj = ccrs.NorthPolarStereo()
         poly = sgeom.Polygon([(0, 0), (-90, 0), (-180, 0), (-270, 0)],