在Python中读取和处理多个netcdf文件

Question

我需要帮助来读取多个netCDF文件，虽然这里有一些例子，但都没有正常工作。我使用的是Python(x,y)版本2.7.5，还有其他一些包：netcdf4 1.0.7-4，matplotlib 1.3.1-4，numpy 1.8，pandas 0.12，basemap 1.0.2……

我有一些在GrADS中习惯做的事情，现在需要在Python中实现。我有一些2米高度的温度数据（每4小时一次的数据，每年来自ECMWF），每个文件包含2米温度数据，X大小为480，Y大小为241，Z大小（层数）为1，T大小（时间）为1460或1464（闰年）。这些文件的名字大致是这样的：t2m.1981.nc，t2m.1982.nc，t2m.1983.nc……等等。

根据这个页面：循环读取netcdf文件并进行计算 - Python或R，这是我现在的进展：

from pylab import *
import netCDF4 as nc
from netCDF4 import *
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np

f = nc.MFDataset('d:/data/ecmwf/t2m.????.nc') # as '????' being the years
t2mtr = f.variables['t2m']

ntimes, ny, nx = shape(t2mtr)
temp2m = zeros((ny,nx),dtype=float64)
print ntimes
for i in xrange(ntimes):
    temp2m += t2mtr[i,:,:] #I'm not sure how to slice this, just wanted to get the 00Z values.
      # is it possible to assign to a new array,...
      #... (for eg.) the average values of  00z for January only from 1981-2000? 

#creating a NetCDF file
nco = nc.Dataset('d:/data/ecmwf/t2m.00zJan.nc','w',clobber=True)
nco.createDimension('x',nx)
nco.createDimension('y',ny)

temp2m_v = nco.createVariable('t2m', 'i4',  ( 'y', 'x'))
temp2m_v.units='Kelvin'
temp2m_v.long_name='2 meter Temperature'
temp2m_v.grid_mapping = 'Lambert_Conformal' # can it be something else or ..
#... eliminated?).This is straight from the solution on that webpage.

lono = nco.createVariable('longitude','f8')
lato = nco.createVariable('latitude','f8')
xo = nco.createVariable('x','f4',('x')) #not sure if this is important
yo = nco.createVariable('y','f4',('y')) #not sure if this is important
lco = nco.createVariable('Lambert_Conformal','i4') #not sure

#copy all the variable attributes from original file
for var in ['longitude','latitude']:
    for att in f.variables[var].ncattrs():
        setattr(nco.variables[var],att,getattr(f.variables[var],att))

# copy variable data for lon,lat,x and y
lono=f.variables['longitude'][:]
lato=f.variables['latitude'][:]
#xo[:]=f.variables['x']
#yo[:]=f.variables['y']

#  write the temp at 2 m data
temp2m_v[:,:]=temp2m

# copy Global attributes from original file
for att in f.ncattrs():
    setattr(nco,att,getattr(f,att))

nco.Conventions='CF-1.6' #not sure what is this.
nco.close()

#attempt to plot the 00zJan mean
file=nc.Dataset('d:/data/ecmwf/t2m.00zJan.nc','r')
t2mtr=file.variables['t2m'][:]
lon=file.variables['longitude'][:]
lat=file.variables['latitude'][:]
clevs=np.arange(0,500.,10.)
map =   Basemap(projection='cyl',llcrnrlat=0.,urcrnrlat=10.,llcrnrlon=97.,urcrnrlon=110.,resolution='i')
x,y=map(*np.meshgrid(lon,lat))
cs = map.contourf(x,y,t2mtr,clevs,extend='both')
map.drawcoastlines()
map.drawcountries()
plt.plot(cs)
plt.show()

第一个问题出现在temp2m += t2mtr[1,:,:]这行。我不确定怎么切片数据，只获取所有文件中的00z（比如说只获取一月份的数据）。

第二个问题是，在运行测试时，cs = map.contourf(x,y,t2mtr,clevs,extend='both')这行出现了错误，提示“形状与z不匹配：发现(1,1)而不是(241,480)”。我知道输出数据可能有错误，可能是记录值时出错，但我搞不清楚是什么问题/在哪里出错。

谢谢你的时间。我希望这不会让人困惑。

Answer 1

所以，t2mtr 是一个三维数组。

ntimes, ny, nx = shape(t2mtr)

这段代码会对第一个维度上的所有值进行求和：

for i in xrange(ntimes):
    temp2m += t2mtr[i,:,:]

有个更好的方法可以做到这一点：

temp2m = np.sum(tm2tr, axis=0)
temp2m = tm2tr.sum(axis=0) # alt

如果你想要平均值，可以用 np.mean 代替 np.sum。

如果你想对一部分时间，比如 jan_times，进行平均，可以用这样的表达式：

jan_avg = np.mean(tm2tr[jan_times,:,:], axis=0)

如果你只想要一个简单的范围，比如前30个时间点，这样做是最简单的。为了简单起见，我假设数据是按天记录的，年份长度是固定的。你可以根据4小时的频率和闰年进行调整。

tm2tr[0:31,:,:]

获取几年的1月数据的简单方法是构建一个索引，比如：

yr_starts = np.arange(0,3)*365 # can adjust for leap years
jan_times = (yr_starts[:,None]+ np.arange(31)).flatten()
# array([  0,   1,   2, ...  29,  30, 365, ..., 756, 757, 758, 759, 760])

另一种选择是重新调整 tm2tr 的形状（但这对闰年不太适用）。

tm2tr.reshape(nyrs, 365, nx, ny)[:,0:31,:,:].mean(axis=1)

你可以用这样的方式测试时间采样：

np.arange(5*365).reshape(5,365)[:,0:31].mean(axis=1)

数据集里没有时间变量吗？你可能可以从中提取出所需的时间索引。我几年前处理过ECMWF的数据，但不太记得细节了。

至于你的 contourf 错误，我建议检查一下三个主要参数的形状：x、y、t2mtr。它们的形状应该是匹配的。我没有使用过 Basemap。