2014-02-19 3 views
1

Я загрузил некоторые данные батиметрии gebco в виде файла netCDF. Я хотел бы построить его с python-basemap. Я пробовал,plot gebco data in python baseemap

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 map

+0

Какова мотивация создания 'X' и' Y' (например, вы делаете это дважды) или для 'Z = zz [:]. Reshape (shape (X))'? Вы также не определили 'nx' или' ny' в предоставленном коде. Что они представляют? –

+0

Простите, это было немного неясно. Код обновлен. zz преобразуется, так как это один вектор столбца. – Dave

+0

Что такое 'x' и' y'? Если они представляют собой массивы долготы и широты соответственно, нет необходимости использовать 'makegrid' с' nx' и 'ny'; просто сделайте 'x, y = m (* np.meshgrid (x, y))'. Возможно, это также устраняет проблему? –

ответ

2

так только для записи, вот ответ, который работает, используя комментарии выше:

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

, который производит этот участок. enter image description here