我已经下载了一些gebco 测深数据作为netCDF 文件。我想用 python-basemap 绘制它。我已经尝试过,
import netCDF4
from mpl_toolkits.basemap import Basemap
# Load data
dataset = netCDF4.Dataset('/home/david/Desktop/GEBCO/gebco_08_-30_45_5_65.nc')
# Extract variables
x = dataset.variables['x_range']
y = dataset.variables['y_range']
spacing = dataset.variables['spacing']
# Data limits
nx = (x[-1]-x[0])/spacing[0] # num pts in x-dir
ny = (y[-1]-y[0])/spacing[1] # num pts in y-dir
# Reshape data
zz = dataset.variables['z']
Z = zz[:].reshape(ny, nx)
# setup basemap.
m = Basemap(llcrnrlon=-30,llcrnrlat=45.0,urcrnrlon=5.0,urcrnrlat=65.0,
resolution='i',projection='stere',lon_0=-15.0,lat_0=55.0)
# Set up grid
lons, lats = m.makegrid(nx, ny)
x, y = m(lons, lats)
m.contourf(x, y, flipud(Z))
m.fillcontinents(color='grey')
m.drawparallels(np.arange(10,70,10), labels=[1,0,0,0])
m.drawmeridians(np.arange(-80, 5, 10), labels=[0,0,0,1])
这给出了下图,显然不正确。问题源于区域的定义方式。对于 basemap 区域,由左下角纬度、经度和右上角纬度、经度定义。但gebco 数据采用沿中心线定义的最大和最小经度/纬度。 任何人都有gebco数据的经验或看到解决方案吗?
谢谢 D
最佳答案
因此,仅供记录,这是使用上面的评论的有效答案:
import netCDF4
from mpl_toolkits.basemap import Basemap
# Load data
dataset = netCDF4.Dataset('/usgs/data1/rsignell/bathy/gebco_08_-30_-45_5_65.nc')
# Extract variables
x = dataset.variables['x_range']
y = dataset.variables['y_range']
spacing = dataset.variables['spacing']
# Compute Lat/Lon
nx = (x[-1]-x[0])/spacing[0] # num pts in x-dir
ny = (y[-1]-y[0])/spacing[1] # num pts in y-dir
lon = np.linspace(x[0],x[-1],nx)
lat = np.linspace(y[0],y[-1],ny)
# Reshape data
zz = dataset.variables['z']
Z = zz[:].reshape(ny, nx)
# setup basemap.
m = Basemap(llcrnrlon=-30,llcrnrlat=45.0,urcrnrlon=5.0,urcrnrlat=65.0,
resolution='i',projection='stere',lon_0=-15.0,lat_0=55.0)
x,y = m(*np.meshgrid(lon,lat))
m.contourf(x, y, flipud(Z));
m.fillcontinents(color='grey');
m.drawparallels(np.arange(10,70,10), labels=[1,0,0,0]);
m.drawmeridians(np.arange(-80, 5, 10), labels=[0,0,0,1]);
产生这个图。
关于python - 在Python basemap 中绘制gebco数据,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/21889079/