在 python 底图中绘制 gebco 数据答案

【问题标题】：plot gebco data in python basemap在 python 底图中绘制 gebco 数据
【发布时间】：2014-02-19 18:21:50
【问题描述】：

我已将一些 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])

这给出了下图，显然不正确。问题源于这些区域的定义方式。对于底图区域，由左下角 lat,lon 和右上角 lat,lon 定义。但是 gebco 数据采用沿中心线定义的最大和最小 lon/lat。有人对 gebco 数据有任何经验或看到解决方案吗？

谢谢 D

【问题讨论】：

创建X 和Y（顺便说一句，你做了两次）或Z=zz[:].reshape(shape(X)) 的动机是什么？您还没有在提供的代码中定义nx 或ny。它们代表什么？
抱歉，有点不清楚。代码已更新。 zz 被重新整形，因为它是一个单列向量。
x 和 y 是什么？如果分别是经度和纬度数组，则makegrid与nx和ny无需使用；就做x,y = m(*np.meshgrid(x, y))。也许这也能解决问题？
x 和 y 包含 lon 和 lat 的最大值和最小值，间距是网格步长。我从它们创建了 lon lat 数组。使用您建议的数组而不是 makegrid 效果很好，谢谢。

标签： python netcdf matplotlib-basemap

【解决方案1】：

所以只是为了记录，这是有效的答案，使用上面的 cmets：

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]);

产生这个情节。

【讨论】：