diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py index 500057ade..148d0036e 100644 --- a/lib/cartopy/crs.py +++ b/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: diff --git a/lib/cartopy/tests/test_polygon.py b/lib/cartopy/tests/test_polygon.py index 0b5ecb25b..9b06573f5 100644 --- a/lib/cartopy/tests/test_polygon.py +++ b/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):