Adaptive line resampling

[greglucas]

Inspired by d3-geo and the way they recursively subdivide their line segments to make them smooth-enough. The difference is that they are always working in spherical geometries, whereas we are working from planar to planar, so our checks and splitting of geometries are slightly different. We already had most of the logic for splitting, so we can re-use some of that and only update the recursive mid-point finding. The previous find a boundary point bisection logic is the same and then we continue to do the recursive midpoint resampling after finding the boundary crossing.

The major improvement here is projection symmetry, so that any projected path A->B is identical output to a path B->A. This has always annoyed me and been on my list of things to tackle for quite a while.

I added an example to try and demonstrate the various thresholds and how the bisection looks with different selections.

AI Disclosure

I used AI to help me in this implementation. The general idea was my own as inspired by d3 noted in issues linked here from before AI times. I have reviewed all code and logic produced and updated comments and logic flow myself so that it made sense to me going through everything.

Testing

This updates a lot of test images which I've kept as a separate commit for ease of comparison and updates. Most of them are minor shifts in the points (making them evenly subdivide the path rather than larger jumps previously), but hardly even noticeable and just barely pushing over the threshold. Others like lib/cartopy/tests/mpl/baseline_images/mpl/test_axes/geoaxes_set_boundary_clipping.png looked worse to me at first glance, but I think I've convinced myself it is actually more correct now. The top previously looked like it was getting filled in with the filled contours, but I think that was an artifact from it being a different point selection than it is now. Now, it is symmetric on the top and bottom, which looks nicer to me and likely more "correct".

test_polygon.py has several tests in it that are somewhat hard to reason about unfortunately. They are lines as polygons, so have zero area, then have a hole that is a single line too. I'm not exactly sure what we want the result to be in those cases. I tried to keep everything similar to before, but some of the areas and points being close to boundary edges threw things off. I did go through the

Previous discussions and related issues

closes #1039

This is related to #498, specifically this comment #498 (comment)
but it is separate and not exactly the topology transform being discussed at the top level._

[greglucas]

Hm... I'm not sure why the mac tests are failing. I'm wondering if this has to do with some floating point precision calculations.

Do we need 2+ holes, or can it sometimes be merged into a single hole and still be OK? I don't have a mac to test this on anymore :(
https://github.com/SciTools/cartopy/actions/runs/23387205255/job/68036057213?pr=2647#step:8:120

[rcomer]

I think that test was from this example #2470 (comment)

So one square hole and another hole around the south pole

[greglucas]

OK, it looks like this is a consequence of the new PROJ 9.7.1, versus PROJ 9.5.1 in most of the runners. I was able to test it locally using conda.

Looking at the test and the old examples, I actually don't think the test was testing what we wanted and we were getting a bit lucky. If I print out all of the areas of the holes in that test I get 8 holes with very small area and one hole with a large area. The small holes are all from the old "seam" points along the antimeridian not meeting at exactly the same spot and producing slight holes.

In the test, we call project_geometry directly, which actually calls .to_geodesic() on the projection and turns it into a spherical interpolation, so the start/end points are the same and there is no resampling done. Also because we are at the pole, we'd need to add more than a few other points to get the desired geodesic, otherwise it would short-circuit back over the pole because that is the shortest path.

cartopy/lib/cartopy/tests/test_polygon.py

Lines 453 to 455 in 05e5c20

	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)

I think we'll need to update the test to use PlateCarree() explicitly, but then I don't think it will be a "hole" in the polar region, but rather two separate polygons split at the seam. There are a few things this did catch though, like I am currently splitting the spherical interpolation points, which doesn't make sense since it is a continuous geometry in that domain. I will have another look at some of this and what is going on in more detail to make sure I'm following all the expectations.