我正在尝试制作一些海冰数据图。数据在 EASE-North 网格中交付,示例文件 (HDF4) 可以在以下位置下载:
ftp://n4ftl01u.ecs.nasa.gov/SAN/OTHR/NISE.004/2013.09.30/
我为 EASE-Grid 创建了一个自定义投影类,它似乎可以正常工作(海岸线与数据对齐得很好)。
当我尝试添加自然地球特征时,它返回一个空的 Matplotlib 图。
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())
上面的脚本工作正常并产生了这个结果:
但是当我取消对 ax.add_feature(land)
的注释时,它失败了,没有任何错误,只返回空图。我是否遗漏了一些明显的东西?
这是 IPython 笔记本: http://nbviewer.ipython.org/6779935
我的 Cartopy 版本是 Christoph Gohlke 网站上的 0.9 版(谢谢!)。
编辑:
尝试保存图形会抛出异常:
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
检查“土地”cartopy.feature
没有发现任何问题,所有多边形都通过了 .isvalid()
并且所有环(ext en int)都为 4 或更多元组。所以输入形状似乎不是问题所在(并且在 PlateCaree() 中工作正常)。
也许某些环(如南半球)在转换为 EASE_North 后会“损坏”?
编辑2:
当我删除内置 NE 功能并加载相同的 shapefile(但剪切低于 40N 的任何内容)时,它会起作用。所以这似乎是某种重投影问题。
for state in shpreader.Reader(r'D:\ne_50m_land_clipped.shp').geometries():
ax.add_geometries([state], cartopy.crs.PlateCarree(),facecolor='#cccccc', edgecolor='#cccccc')
最佳答案
我会说这是一个错误。我猜 add_feature
更新了 matplotlib viewLim,结果是图片放大到一个很小的区域(除非你缩小很多,否则它会显示为白色)。
从我的头脑来看,我认为底层行为在 matplotlib 中得到了改进,但 cartopy 尚未使用新的 viewLim 计算。同时,我建议您手动设置 map 的范围:
ax.set_extent(范围,transform=EASE_North())
HTH
关于matplotlib - 在自定义投影上绘制自然地球特征,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/19119778/