我在读取多个 netCDF 文件时需要帮助,尽管这里的示例很少,但它们都无法正常工作。 我正在使用 Python(x,y) vers 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 米温度数据,Xsize=480,Ysize=241, Zsize(level)=1, Tsize(time) = 1460 or 1464 for leap years. 这些是我的文件名看起来很像:t2m.1981.nc、t2m.1982.nc、t2m.1983.nc ...等。
基于此页面: ( Loop through netcdf files and run calculations - Python or 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')
出现错误,提示“shape does not match that of z : 发现 (1,1) 而不是 (241,480)”。我知道输出数据可能存在一些错误,这是由于记录值时出错,但我不知道是什么/在哪里。
感谢您的宝贵时间。我希望这不会造成混淆。
最佳答案
所以 t2mtr
是一个 3d 数组
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,:,:]
获取几年一月数据的一种简单方法是构建如下索引:
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])
另一种选择是 reshape 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
错误,我会检查 3 个主要参数的形状:x
,y
,t2mtr
。他们应该匹配。我没有使用过 Basemap
。
关于python - 在python中读取和操作多个netcdf文件,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/22287039/