Repository: SciTools/cartopy
Pull request #2521 (2025-04-23 16:23:56)

FIX: make hole multipolygon valid before removing from projected polygon
<!--

Thanks for contributing to cartopy!
Please use this template as a guide to streamline the pull request you are about to make.

Remember: it is significantly easier to merge small pull requests. Consider whether this pull
request could be broken into smaller parts before submitting.

-->



## Rationale

<!-- Please provide detail as to *why* you are making this change. -->
When a polygon with multiple holes is projected, the projected holes may become so close that they can't form a valid MultiPolygon as they [touch along a line](https://shapely.readthedocs.io/en/stable/manual.html#MultiPolygon).  In this case, we can still create the holes one by one.  Taking the first example from https://github.com/SciTools/cartopy/issues/2176#issuecomment-2772601094, and adding coastlines for reference, 

```python
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt

crs_proj = ccrs.LambertAzimuthalEqualArea(central_latitude=-90)
fig, ax = plt.subplots(1, 1, subplot_kw={"projection": crs_proj})
ax.add_feature(cfeature.OCEAN)
ax.coastlines()
plt.show()
```
I now get
![image](https://github.com/user-attachments/assets/b5635630-c98c-4158-ab62-1f5237e15d37)

Setting the central latitude to -30 still results in a failure, but this is very sensitive to the exact number and the user will be able to work around it by setting central latitude to e.g. -29.9999:

![image](https://github.com/user-attachments/assets/a7cfab2e-ea94-493d-9354-b5edbcfa0c29)


## Implications

<!-- If applicable, to the best of your knowledge, what are the implications of this change? -->
Making the holes one by one is slow compared with using the MultiPolygon:  I tried removing the try-except and always removing them in the loop.  This benchmark took 18 times longer:

```python
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature


def ocean_fig():
    fig = plt.figure()
    ax = fig.add_subplot(111, projection=ccrs.NorthPolarStereo())

    ax.add_feature(cfeature.OCEAN.with_scale('10m'))
    
    return fig


import timeit

setup = """
from __main__ import ocean_fig
fig = ocean_fig()
"""

print(timeit.timeit('fig.draw_without_rendering()', setup=setup, number=5))  
```

<!--
## Checklist

 * If this is a new feature, please provide an example of its use in the description. We may want to make a
   follow-on pull request to put the example in the gallery!

 * Ensure there is a suitable item in the cartopy test suite for the change you are proposing.

 * Note that you will automatically be asked to sign the
   [Contributor Licence Agreement](https://cla-assistant.io/SciTools/)
   (CLA), if you have not already done so.

-->


rcomer: CI failures are in `test_wfs_france` so not my fault!
greglucas: > CI failures are in test_wfs_france so not my fault!

I noticed that on main too. We should probably xfail that for now.

Patch:
diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py
index 07fce8e0..f58dc2ad 100644
--- a/lib/cartopy/crs.py
+++ b/lib/cartopy/crs.py
@@ -1251,7 +1251,8 @@ class Projection(CRS, metaclass=ABCMeta):
             box = sgeom.box(x3, y3, x4, y4, ccw=is_ccw)
 
             # Invert the polygons
-            polygon = box.difference(sgeom.MultiPolygon(interior_polys))
+            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)


Test patch:
diff --git a/lib/cartopy/tests/test_polygon.py b/lib/cartopy/tests/test_polygon.py
index 22a3ee57..c9fb435a 100644
--- a/lib/cartopy/tests/test_polygon.py
+++ b/lib/cartopy/tests/test_polygon.py
@@ -457,6 +457,14 @@ class TestHoles(PolygonTests):
         assert len(multi_polygon.geoms) == 1
         assert len(multi_polygon.geoms[0].interiors) >= 2
 
+    def test_inverted_poly_merged_holes(self):
+        proj = ccrs.LambertAzimuthalEqualArea(central_latitude=-90)
+        pc = ccrs.PlateCarree()
+        poly = sgeom.Polygon([(-180, -80), (180, -80), (180, 90), (-180, 90)],
+                             [[(-50, 60), (-50, 80), (0, 80), (0, 60)],
+                              [(-50, 81), (-50, 85), (0, 85), (0, 81)]])
+        # Smoke test that nearby holes do not cause side location conflict
+        proj.project_geometry(poly, pc)
 
     def test_inverted_poly_clipped_hole(self):
         proj = ccrs.NorthPolarStereo()