FIX: create a single inverted polygon when no exteriors found

[rcomer]

Rationale

When _rings_to_multipolygon finds interior rings without associated exteriors, for each of these interiors it creates a new polygon as the difference between ~the whole map and that interior. Intuitively, if there are multiple interiors, then all of those whole map polygons will overlay each other and you won't see the holes. Instead it should create a single ~whole map polygon which uses all of those left over interiors.

I also switched out use of .buffer(0) for shapely.make_valid as that seems like the more official way to do that.

Fixes #1149, fixes #1674, fixes #2160 (all ocean feature examples) and fixes #1449 (contourf example)
Also fixes the original case from #2030 (but not the second example so not linking that issue)
Also fixes the contourf examples at #1076 (comment) and #1076 (comment)

I turned #1149 into an image test. Really there should be something more specific in test_polygon.py, but I'm at a loss how to create an artificial example for this.

Implications

[rcomer]

Note that MultiPolygon will work fine with an empty list.

[rcomer]

We do get a GeometryCollection for the #1149 case.

[rcomer]

I note that _rings_to_multipolygon is called within _project_polygon. If you start with a MultiPolygon it will just loop through the polygons and call _project_polygon on each, so we might plausibly create multiple overlapping whole map polygons that way. So (plausibly) we might fix more cases by applying this logic at the MultiPolygon level. I propose to get this fix in first and dig into the MultiPolygon question later though.

[rcomer]

I appear to be reinventing part of #1739...

[greglucas]

Really there should be something more specific in test_polygon.py, but I'm at a loss how to create an artificial example for this.

I agree. Simplifying these down would really help us debug future issues too because I feel like it is confusing what the root cause of the issue is sometimes. A few thoughts in case this sparks anything.

• All geometries causing issues are coming from the default PlateCarree, so we should be able to create something in degrees.
• Most of the issues seem to be around Lambert* transforms, so going to an r/theta coordinate frame.
• Are there situations where the bottom left corner of a box gets transformed into the upper right corner of a box when you are right around the central longitude / maximum latitude? (possibly a rhombus/diamond shape where one point is on the opposite side of a boundary)

I appear to be reinventing part of #1739...

That has a lot going on. Maybe this is a small piece of that and we can keep doing incremental updates?

[rcomer]

OK here's the simplest I have so far by studying #1149 (comment)

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import shapely.geometry as sgeom

sudo_ocean = sgeom.Polygon([(-180, -80), (180, -80), (180, 90), (-180, 90)],
                           [[(-50, -50), (-50, 0), (0, 0), (0, -50)]])

fig = plt.figure(figsize=(6, 8))

ax1 = fig.add_subplot(2, 1, 1, projection=ccrs.PlateCarree())
proj = ccrs.LambertAzimuthalEqualArea(central_latitude=45, central_longitude=-100)
ax2 = fig.add_subplot(2, 1, 2, projection=proj)

for ax in [ax1, ax2]:
    ax.add_geometries(sudo_ocean, crs=ccrs.PlateCarree())
    ax.set_global()
    
plt.show()

main:

branch:

[rcomer]

Should I remove the image test?

[rcomer]

We actually get 36 holes. I assume most of them are responsible for the seam down the dateline visible in #2470 (comment). I didn't assert the exact number because I've no idea how stable that would be.

Against main, this test fails at the first assert since we have two polygons in our multipolygon.

[greglucas]

I would remove the image test since it is downloading the natural earth ocean for every runner, or replace it with this example as an image comparison test that doesn't require external downloading.

The 36 holes is interesting and worth investigating IMO. For your example here, it looks like there are many small loops around the boundary points that are created. (your image has a small white curve/line in it protruding from the south pole area which are these linear rings I think)

import matplotlib.pyplot as plt

ax = plt.figure().add_subplot()

interior = multi_polygon.geoms[0].interiors
for i in range(len(interior)):
    ax.plot(*interior[i].xy)

plt.show()

I wonder if these extra loops/rings are what is slowing down some of our projecting of geometries even 🤔

I'm also happy to defer these investigations to future PRs because this is already a pretty substantial improvement I think.