FIX: handle exterior ring - boundary attachment inversion

Some rings can get inverted when being transformed and the exterior
should actually be classified as an interior hole. Catch that case
for when an interior area is larger than the exterior area and
invert the polygons.

Author: greglucas

lib/cartopy/crs.py

@@ -917,7 +917,10 @@ def _project_linear_ring(self, linear_ring, src_crs):
                                                     line.coords[-1],
                                                     atol=threshold):
                 result_geometry = sgeom.LinearRing(line.coords[:-1])
-                rings.append(result_geometry)
+                # Discard near-zero area rings; they are boundary artefacts
+                # produced when a source vertex lies just outside the domain.
+                if sgeom.Polygon(result_geometry).area >= threshold ** 2:
+                    rings.append(result_geometry)
             else:
                 line_strings.append(line)
         # If we found any rings, then we should re-create the multi-line str.
@@ -1190,6 +1193,26 @@ def _rings_to_multi_polygon(self, rings, is_ccw):
             else:
                 exterior_rings.append(ring)
 
+        # Detect "inverted" exterior rings: when a projected exterior ring is
+        # smaller than an interior ring, the source polygon spans a region
+        # that wraps around the projection (e.g. a global polygon projected
+        # to a bounded CRS). Such a ring represents the polar/antipodal
+        # boundary that should be a hole, not an exterior boundary. Treat
+        # it as an interior ring so the full projection domain is used as the
+        # implicit exterior via the box-difference path below.
+        # Only apply when all exterior rings are smaller than an interior ring,
+        # which signals the whole polygon is "inverted" rather than merely
+        # having one small piece from an antimeridian split.
+        if interior_rings and exterior_rings:
+            interior_areas = [sgeom.Polygon(r).area for r in interior_rings]
+            max_interior_area = max(interior_areas)
+            reclassified = [r for r in exterior_rings
+                            if sgeom.Polygon(r).area < max_interior_area]
+            if len(reclassified) == len(exterior_rings):
+                for r in reclassified:
+                    exterior_rings.remove(r)
+                    interior_rings.append(r)
+
         polygon_bits = []
 
         # Turn all the exterior rings into polygon definitions,

lib/cartopy/tests/test_polygon.py

@@ -453,10 +453,15 @@ def test_inverted_poly_multi_hole(self):
             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
+        pc = ccrs.PlateCarree()
+        multi_polygon = proj.project_geometry(poly, pc)
         assert len(multi_polygon.geoms) == 1
         assert len(multi_polygon.geoms[0].interiors) >= 2
+        # A point below the -80 line shouldn't be covered up by the
+        # projected polygon
+        pt = proj.project_geometry(sgeom.Point(0, -85), pc)
+        assert not multi_polygon.contains(pt)
+
 
     def test_inverted_poly_merged_holes(self):
         proj = ccrs.LambertAzimuthalEqualArea(central_latitude=-90)