在自定义投影上绘制自然地球特征 [英] Plotting Natural Earth features on a custom projection
问题描述
我正在尝试绘制一些海冰数据图.数据以EASE-North网格形式提供，可以从以下位置下载示例文件(HDF4):
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/
我为EASE-Grid创建了一个自定义投影类，它似乎正在起作用(海岸线与数据完全吻合).
I created a custom projection class for the EASE-Grid, it seems to be working (the coastlines align well with the data).
当我尝试添加自然地球特征时，它返回一个空的 Matplotlib 图.
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:
但是当我取消注释 ax.add_feature(land) 时，它失败了，没有任何错误，只返回空图.我是否遗漏了一些明显的东西?
But when i uncomment the ax.add_feature(land) it fails without any error, only returning the empty figure. Am i missing something obvious?
这是IPython Notebook:http://nbviewer.ipython.org/6779935
Here is the IPython Notebook: http://nbviewer.ipython.org/6779935
我的 Cartopy 版本是 Christoph Gohlke 网站上的 0.9 版(谢谢！).
My Cartopy build is version 0.9 from Christoph Gohlke's website (thanks!).
尝试保存图形确实会引发异常:
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
检查土地"cartopy.feature 没有发现任何问题，所有多边形都通过了 .isvalid() 并且所有环 (ext en int) 都是 4 或更多元组.因此，输入形状似乎不是问题所在(并且在PlateCaree()中工作正常).
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()).
也许某些环(如南半球)在转换为 EASE_North 后会损坏"?
Maybe some rings (like on the southern hemisphere) get 'corrupt' after transforming to EASE_North?
当我删除内置NE功能并加载相同的shapefile(但剪辑后的内容小于40N)时，它会起作用.所以这似乎是某种重新投影问题.
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')
推荐答案
我会说这是一个错误.我猜想 add_feature 更新了matplotlib viewLim，结果是图片放大到一个很小的区域(除非您将其缩小很多，否则它将显示为白色).
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).
从头开始，我认为matplotlib的基本行为已得到改善，但是cartopy尚未使用新的viewLim计算.同时，我建议您使用以下方法手动设置地图的范围:
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
这篇关于在自定义投影上绘制自然地球特征的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
