Treatment of interior rings from transform between feature and projection domains

[kdpenner]

Rationale

Fixes #1685, #1613, #1149, #1076, #1723, #955, #1367, #1617, and a bunch of other unfiled issues. NB a few of them were incorrectly marked as duplicates. Methodology applicable to #1362. Potentially solves #323 and #1449?

Some of the test suite failures are expected.

This PR is a work in progress. Debugging junk is everywhere and more work is needed. I'm submitting now as a heads-up; the changes will need extensive testing. In particular I've done no debugging of the sampling of non-Plate Carree feature CRSs and have written no unit tests.

A summary of the discoveries this PR solves:

1. points that project to infinity are improperly included (-> changes to FeatureArtist)
2. _project_linear_rings sometimes creates small interior rings. _rings_to_multi_polygon inverts the rings, which leads to polygons covering almost the entire domain (-> changes to _rings_to_multi_polygon)
3. OTOH _project_linear_rings sometimes creates exterior rings encircling the entire domain or almost the entire domain. If two of these rings are created the interior ring removal logic leads to polygons covering almost the entire domain (-> changes to _rings_to_multi_polygon)
4. correcting 2) or 3) invalidates the MultiPolygon creation process (-> changes to _project_polygon and _project_multipolygon)
5. _project_linear_rings sometimes flips the orientation of the rings. i.e. land polygons are filled when ocean polygons are passed

Implications

Performance is reduced, severely for large polygons. If invalid in their native CRS they have to be buffered for a difference, then buffered again in the projection CRS. But the following plot is successful.

[kdpenner]

Note to self just discovered inversion of coastlines in middle plots of bottom 2 rows

Edit: they're weren't inverted---they weren't added to the dictionary

[dopplershift]

@kdpenner I don't have time right at this moment to understand everything that's going on here, but wanted to say thanks for putting this in. Clearly this goes a long way towards addressing a lot of problematic corner cases we have.

Do you have a benchmark you can post that led you to the "performance is reduced" comment?

Regarding the unit tests, in the absence of some kind of systematic approach we can take here, I think all of the linked issues represent a good source of tests.

[kdpenner]

@dopplershift you're welcome. I hope to have a final version in a month or so. (Depends on my other work.) Will present a benchmark then.

[kdpenner]

Ready for review. I'll have an update with benchmarks in a few days.

[kdpenner]

Also: much of this feels like a treatment of symptoms and not the illness. My guess is that there are flaws in _attach_lines_to_boundary and _project_linear_ring and maybe project_linear, but in those parts of the code there are variables named thing; I don't want to go that far down the call sequence.

[QuLogic]

This appears to need a rebase.

[kdpenner]

Done

[kdpenner]

A simple benchmark script:

import timeit
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt


def tmerc():
    fig = plt.figure()
    proj = ccrs.TransverseMercator(central_longitude=18.14159, approx=False)
    ax = fig.add_subplot(1, 1, 1, projection=proj)
    ax.add_feature(cfeature.OCEAN.with_scale("10m"))
    plt.savefig("1.png", bbox_inches="tight")
    plt.close()


print(timeit.timeit(tmerc, number=1))


oceans = list(cfeature.OCEAN.with_scale("10m").geometries())[0]


def buffer_oceans():
    oceans.buffer(0)


print(timeit.timeit(buffer_oceans, number=1))

runs in 24s with cartopy 0.19. All of the time is in the first run so I didn't use number > 1.

With this PR the script runs in 84s.

Replacing the ocean feature with ax.coastlines():

v0.19: 0.25s
   PR: 0.37s

[QuLogic]

runs in 24s with cartopy 0.19. All of the time is in the first run so I didn't use number > 1.

With this PR the script runs in 84s.

3.5x slower? That seems a pretty hefty penalty. I wonder if there is some way to speed that up.

[kdpenner]

Probably. There's a double for loop, over the exterior and interior rings. There are expensive unary unions. But the major reason for the time increment is a buffer of invalid polygons in their native CRS, and I see no way around it. shapely's make_valid, introduced in 1.8, may be faster than buffer?

I modified the benchmark script in my previous comment; here are notional numbers.

version	action	runtime
0.19	add entire ocean feature	24s
PR	add entire ocean feature	84s
...	buffer entire ocean feature	72s
PR	add entire ocean feature after buffering	10s
0.19	add Caspian Sea (valid polygon)	0.24s
PR	add Caspian Sea (valid polygon)	0.29s
0.19	add coastlines	0.24s
PR	add coastlines	0.32s

Edit: add a line in the table

[CLAassistant]

Thank you for your submission! We really appreciate it. Like many open source projects, we ask that you sign our Contributor License Agreement before we can accept your contribution.
_{You have signed the CLA already but the status is still pending? Let us recheck it.}