SciTools / cartopy

Cartopy - a cartographic python library with matplotlib support
https://scitools.org.uk/cartopy/docs/latest
BSD 3-Clause "New" or "Revised" License
1.41k stars 359 forks source link

Broken ocean polygon reprojection with LambertAzimuthalEqualArea #1149

Open QuLogic opened 5 years ago

QuLogic commented 5 years ago

May I add LambertAzimuthalEqualArea to the list or reproducers with current master + matplotlib2

proj = cartopy.crs.LambertAzimuthalEqualArea(central_latitude=45,
                                             central_longitude=-100)
ax = plt.axes(projection=proj)
ax.add_feature(cartopy.feature.OCEAN)
ax.coastlines()
ax.set_extent([-140, -70, 20, 60])

test

Originally posted by @akrherz in https://github.com/SciTools/cartopy/issues/803#issuecomment-333909742

QuLogic commented 5 years ago

Seems to also affect AzimuthalEquidistant (which is somewhat similar.)

QuLogic commented 5 years ago

If we plot the points of the OCEAN polygons:

for i, mp in enumerate(cartopy.feature.OCEAN.geometries()):
    p = mp[0]  # There's only one polygon.
    x, y = zip(*list(p.exterior.coords))
    ax.scatter(x, y, marker='o', facecolor='none', edgecolor=f'C{i * 2}',
               transform=cartopy.feature.OCEAN.crs,
               label=f'{i} exterior')
    for j, interior in enumerate(p.interiors):
        x, y = zip(*list(interior.coords))
        ax.scatter(x, y, marker='x', color=f'C{i * 2 + 1}',
                   transform=cartopy.feature.OCEAN.crs,
                   label=f'{i} interior' if j == 0 else None)

figure_1 then we see that there's one small 'ocean' in Eurasia and one big ocean whose exterior is Eurasia+Africa and Antarctica, with the remaining continents being the excluded interiors. This is a bit more clear on a PlateCarree(central_longitude=30): figure_1-1 It's a bit weird to me that both Eurasia and Antarctica are both part of the exterior polygon, but I think something weird is going on there, either tranform-wise, or attaching to the boundary.

QuLogic commented 5 years ago

So I think the problem here is re-projecting the weird exterior produces an invalid LinearRing. Using this simplified topology: exterior-pc the re-projection looks like: exterior-laea and the bit extending from the Eastern part of Asia down to Antarctica crosses back over Asia. I think this is an invalid LinearRing to shapely. It's a bit weird to me that that whole line does not seem to be problematic on say, shifted PlateCarree, though.

QuLogic commented 5 years ago

This is essentially the same as #146.