Fix inside-out polygon artifact in ObliqueMercator projection

lib/cartopy/crs.py

@@ -1192,10 +1192,38 @@ def _rings_to_multi_polygon(self, rings, is_ccw):
 
         polygon_bits = []
 
+        # Precompute the projection boundary polygon for inversion detection.
+        boundary_poly = self.domain
+
         # Turn all the exterior rings into polygon definitions,
         # "slurping up" any interior rings they contain.
         for exterior_ring in exterior_rings:
             polygon = sgeom.Polygon(exterior_ring)
+
+            # When a polygon crosses the projection cut line,
+            # _attach_lines_to_boundary can produce a ring that traces the
+            # *complement* of the intended shape while still carrying the
+            # correct CCW/CW winding flag (so it is classified as an exterior
+            # ring rather than an interior one).  The tell-tale sign is an
+            # area that exceeds half the projection domain: a genuine land or
+            # ocean polygon can never be that large.  Invert such rings back
+            # to the correct shape before continuing.
+            if polygon.area > 0.5 * boundary_poly.area:
+                corrected = boundary_poly.difference(
+                    shapely.make_valid(polygon))
+                if not corrected.is_empty:
+                    if isinstance(corrected, sgeom.MultiPolygon):
+                        for p in corrected.geoms:
+                            if not p.is_empty:
+                                polygon_bits.append(
+                                    (p.exterior.coords,
+                                     [i.coords for i in p.interiors]))
+                    elif isinstance(corrected, sgeom.Polygon):
+                        polygon_bits.append(
+                            (corrected.exterior.coords,
+                             [i.coords for i in corrected.interiors]))
+                continue
+
             prep_polygon = prep(polygon)
             holes = []
             for interior_ring in interior_rings[:]:
@@ -1213,7 +1241,6 @@ def _rings_to_multi_polygon(self, rings, is_ccw):
         # Any left over "interior" rings need "inverting" with respect
         # to the boundary.
         if interior_rings:
-            boundary_poly = self.domain
             x3, y3, x4, y4 = boundary_poly.bounds
             bx = (x4 - x3) * 0.1
             by = (y4 - y3) * 0.1
@@ -1257,11 +1284,8 @@ def _rings_to_multi_polygon(self, rings, is_ccw):
             box = sgeom.box(x3, y3, x4, y4, ccw=is_ccw)
 
             if interior_polys:
-                # Invert any valid interior polygons
                 multi_poly = shapely.make_valid(sgeom.MultiPolygon(interior_polys))
                 polygon = box.difference(multi_poly)
-
-                # Intersect the inverted polygon with the boundary
                 polygon = boundary_poly.intersection(polygon)
 
                 if not polygon.is_empty:

lib/cartopy/tests/test_polygon.py

@@ -281,6 +281,39 @@ def test_project_degenerate_poly(self):
         polygons = target.project_geometry(polygon, source)
         assert isinstance(polygons, sgeom.MultiPolygon)
 
+    @pytest.mark.natural_earth
+    def test_oblique_mercator_no_inverted_land(self):
+        # Regression test for inside-out polygon artifact in ObliqueMercator.
+        # When a land polygon crosses the projection cut line,
+        # _rings_to_multi_polygon previously misidentified the winding order
+        # and returned the complement of the intended shape — a polygon
+        # covering the entire projection domain minus a small hole.
+        # These inverted polygons caused ocean areas to render as land.
+        # See https://github.com/SciTools/cartopy/issues/XXXX
+        import cartopy.feature as cfeature
+
+        projection = ccrs.ObliqueMercator(
+            central_longitude=-161.07,
+            central_latitude=54.55,
+            azimuth=90.0,
+            scale_factor=1,
+        )
+        boundary_area = sgeom.Polygon(projection.boundary).area
+        src_crs = ccrs.PlateCarree()
+
+        # The bug is only triggered with 50m (or finer) Natural Earth data.
+        # At small map extents, AdaptiveScaler selects 50m data, but the
+        # 110m dataset does not contain the affected polygons.
+        for geom in cfeature.LAND.with_scale('50m').geometries():
+            projected = projection.project_geometry(geom, src_crs)
+            if not projected.is_empty:
+                assert projected.area <= 0.5 * boundary_area, (
+                    f'project_geometry returned an inverted polygon '
+                    f'(area {projected.area:.3e} > 0.5 * boundary '
+                    f'{0.5 * boundary_area:.3e}). This is the inside-out '
+                    f'polygon artifact where a land polygon crossing the '
+                    f'projection cut line is rendered as its complement.')
+
 
 class TestQuality:
     def setup_class(self):