Oceans (10m resolution) fill the whole plot when using TransverseMercator with central_longitude=15.0

[dbratell]

Description

When plotting oceans (resolution "10m") with projection TransverseMercator(central_longitude=15.0), the whole plot turns blue.

I have understood that ocean polygons are complicated and upgrading to 0.18 made resolution "50m" work, but "10m" still causes the whole plot to be blue.

I've tested various values for central_longitude. There is no bug for values between 10.0 and 12.0. Everything is blue between 12.5 and 19.0, except for 18.0 that for some reason seems to work. I've also tested approx=False and approx=True and that makes no observable difference.

TransverseMercator(central_longitude=15.0) (with some additional parameters) is the official projection for Sweden which makes that particular value interesting.

Code to reproduce

A minimal testcase:

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

MAP_PROJECTION = ccrs.TransverseMercator(central_longitude=15.0)
ax = plt.axes(projection=MAP_PROJECTION)
ax.add_feature(cfeature.OCEAN.with_scale("10m")) # "50m" works
plt.show()

Expected,

The Baltic Sea, and the North Sea

Actually

Blue

### Operating system

Windows 10

Cartopy version

cartopy.version is '0.18.0'
Installed with anaconda/conda which calls the package 0.18.0

Python version

(base) c:>python --version
Python 3.7.8

[dopplershift]

Looks similar to behavior seen in #1578 .

I appreciate the detailed testing here. Whoever debugs this (maybe me eventually), will probably start needing to go through polygon by polygon to see what's failing and why.

Out of curiosity, do you have fiona installed? In #1578, the error was solved/worked-around by installing fiona and relying on that to read in the shapefile data rather than pyshp.

[dbratell]

I've not used fiona before so maybe I got this wrong, but it looks to me like I got the exact same result:

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import fiona
import matplotlib.pyplot as plt
import shapely

with fiona.open("ne_10m_ocean/ne_10m_ocean.shp") as oceans:
    geometries = [shapely.geometry.shape(x["geometry"]) for x in oceans]

MAP_PROJECTION = ccrs.TransverseMercator(central_longitude=15.0, approx=False)
ax = plt.subplot(1, 1, 1, projection=MAP_PROJECTION)
ax.add_feature(cfeature.ShapelyFeature(geometries, ccrs.PlateCarree()))
plt.show()

By the way, when I said that I expected the Baltic sea, that was when I had this line included
ax.set_extent((2, 50, 54, 66)) # Scandinavia and the Baltic Sea

I had removed it from the minimal testcase so with the minimal testcase you would expect to see most of the world instead.

Anyway, looking like fiona doesn't help, unless I did something wrong above.

[dbratell]

I don't know if it's of any help, but when trying all possible non-very-distorting projections to find an alternative projection to use, I also reproduced the problem with ccrs.AzimuthalEquidistant(central_longitude=15, central_latitude=62) and ccrs.LambertAzimuthalEqualArea(central_longitude=15, central_latitude=62). Many others were fine.

[kdpenner]

This is a problem for more than the ocean shapefiles.

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

latlon_crs_plt = ccrs.PlateCarree()
tranmerc_crs_plt = ccrs.TransverseMercator(central_longitude=0., approx=False)

fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot(1, 1, 1, projection=tranmerc_crs_plt)

ax.coastlines()

plt.show()

produces:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

latlon_crs_plt = ccrs.PlateCarree()
tranmerc_crs_plt = ccrs.TransverseMercator(central_longitude=15., approx=False)

fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot(1, 1, 1, projection=tranmerc_crs_plt)

ax.coastlines()

plt.show()

produces:

[kdpenner]

tl;dr

There are rings and holes in lat-lon space that transverse Mercator projects to infinity. Geometry coordinates which fall in the rings and holes are projected to infinity. I suspect cartopy improperly handles infinite coordinates when constructing geometries in the projected space.

A few things going on here

1. there are two ocean polygons at 10m scale in the natural earth dataset. One is invalid. buffer(0) black magic makes it valid.
2. the leftmost column of the figure shows one set of rings and holes that project to infinity.
3. cartopy checks only that the ocean polygons intersect the global extent
   

   cartopy/lib/cartopy/feature/__init__.py

   Lines 103 to 109 in efd15f1

   	if extent is not None:
   	extent_geom = sgeom.box(extent[0], extent[2],
   	extent[1], extent[3])
   	return (geom for geom in self.geometries() if
   	geom is not None and extent_geom.intersects(geom))
   	else:
   	return self.geometries()

   
   For this check the ocean polygons can be invalid and any of the coordinates can project to infinity.
4. ideally you'd construct the lat-lon space that projects to infinity and subtract it from the ocean polygons. Difference requires valid polygons, hence the need to buffer(0).
5. I couldn't find an analytic formulation for the rings and holes, so I fudged them using the regular grid in the leftmost plots.
6. the results of cleaning the ocean polygons, subtracting the regions going to Infinity and Beyond, and projecting the difference are shown in the rightmost plots.
7. Sweden gets a waterworld because exterior polygon coordinates corresponding to Indonesia fall in the rings and holes for longitude_natural_origin = 15.
8. I fixed the coastlines issue by following more or less the same procedure. Lines are simpler than polygons so I used a different implementation.

[greglucas]

@kdpenner, PRs are certainly welcome!

I'm not sure if this helps, but this was a pretty good start of fixing the regions/holes if you want to dig into that code as well.
#1479 (comment)

[kdpenner]

@greglucas Sure, I'll look at this. I don't know how generalizable my hack is---I've tuned some of the parameters---though it does look similar to the one in the PR comment.