import gdal
import cartopy
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)
# example data from:
# ftp://n4ftl01u.ecs.nasa.gov/SAN/OTHR/NISE.004/2013.09.30/
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)
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())
<matplotlib.image.AxesImage at 0x3f9f7b0>