Plotting Natural Earth features on a custom projection

Question

I am trying to make some plots of sea ice data. The data is delivered in the EASE-North grid, an example file (HDF4) can be downloaded at:

ftp://n4ftl01u.ecs.nasa.gov/SAN/OTHR/NISE.004/2013.09.30/

I created a custom projection class for the EASE-Grid, it seems to be working (the coastlines align well with the data).

When i try to add a Natural Earth feature, it returns an empty Matplotlib figure.

import gdal
import cartopy

# projection class
class EASE_North(cartopy.crs.Projection):
    
    def __init__(self):
        
        # see: http://www.spatialreference.org/ref/epsg/3408/
        proj4_params = {'proj': 'laea',
            'lat_0': 90.,
            'lon_0': 0,
            'x_0': 0,
            'y_0': 0,
            'a': 6371228,
            'b': 6371228,
            'units': 'm',
            'no_defs': ''}
        
        super(EASE_North, self).__init__(proj4_params)
        
    @property
    def boundary(self):
        coords = ((self.x_limits[0], self.y_limits[0]),(self.x_limits[1], self.y_limits[0]),
                  (self.x_limits[1], self.y_limits[1]),(self.x_limits[0], self.y_limits[1]),
                  (self.x_limits[0], self.y_limits[0]))
        
        return cartopy.crs.sgeom.Polygon(coords).exterior
        
    @property
    def threshold(self):
        return 1e5
    
    @property
    def x_limits(self):
        return (-9000000, 9000000)

    @property
    def y_limits(self):
        return (-9000000, 9000000)


# read the data
ds = gdal.Open('D:/NISE_SSMISF17_20130930.HDFEOS')

# this loads the layers for both hemispheres
data = np.array([gdal.Open(name, gdal.GA_ReadOnly).ReadAsArray()
                 for name, descr in ds.GetSubDatasets() if 'Extent' in name])

ds = None

# mask anything other then sea ice
sea_ice_concentration = np.ma.masked_where((data < 1) | (data > 100), data, 0)

# plot
lim = 3000000

fig, ax = plt.subplots(figsize=(8,8),subplot_kw={'projection': EASE_North(), 'xlim': [-lim,lim], 'ylim': [-lim,lim]})

land = cartopy.feature.NaturalEarthFeature(
            category='physical',
            name='land',
            scale='50m',
            facecolor='#dddddd',
            edgecolor='none')

#ax.add_feature(land)
ax.coastlines()

# from the metadata in the HDF
extent = [-9036842.762500, 9036842.762500, -9036842.762500, 9036842.762500]

ax.imshow(sea_ice_concentration[0,:,:], cmap=plt.cm.Blues, vmin=1,vmax=100,
          interpolation='none', origin='upper', extent=extent, transform=EASE_North())

The script above works fine and produces this result:

But when i uncomment the ax.add_feature(land) it fails without any error, only returning the empty figure. Am i missing something obvious?

Here is the IPython Notebook: http://nbviewer.ipython.org/6779935

My Cartopy build is version 0.9 from Christoph Gohlke's website (thanks!).

edit:

Trying to save the figure does throw an exception:

fig.savefig(r'D:\test.png')

C:\Python27\Lib\site-packages\shapely\speedups\_speedups.pyd in shapely.speedups._speedups.geos_linearring_from_py (shapely/speedups/_speedups.c:2270)()

ValueError: A LinearRing must have at least 3 coordinate tuples

Examining the 'land' cartopy.feature reveals no issues, all polygons pass the .isvalid() and all rings (ext en int) are of 4 or more tuples. So the input shape doesnt seem to be the problem (and works fine in PlateCaree()).

Maybe some rings (like on the southern hemisphere) get 'corrupt' after transforming to EASE_North?

edit2:

When i remove the build-in NE features and load the same shapefile (but with anything below 40N clipped) it works. So it seems like some sort of reprojection issue.

for state in shpreader.Reader(r'D:\ne_50m_land_clipped.shp').geometries():
    ax.add_geometries([state], cartopy.crs.PlateCarree(),facecolor='#cccccc', edgecolor='#cccccc')

Answer 1

I'd have said that this was a bug. I'm guessing add_feature updates the matplotlib viewLim and the result is that the picture zooms in to a tiny area (which appears white unless you zoom out a lot).

From the top of my head, I think the underlying behaviour has been improved in matplotlib, but cartopy is not yet making use of the new viewLim calculation. In the meantime I'd suggest setting the extents of your map manually with:

ax.set_extent(extent, transform=EASE_North())

HTH