Чтение и управление несколькими файлами netcdf в python

Мне нужна помощь в чтении нескольких файлов 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 метра, с Xsize = 480, Ysize = 241, Zsize (уровень) = 1, Tsize (время) = 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)». Я знаю некоторую ошибку, вероятно, в выходных данных из-за ошибки при записи значений, но я не могу понять, что/где.

Спасибо за ваше время. Надеюсь, это не сбивает с толку.


person Fir Nor    schedule 09.03.2014    source источник
comment
Что вы подразумеваете под данными 00Z? Я не понимаю, в каком формате ваш набор данных: 3D, с формой = (время, x, y), что я понимаю. Где/что такое Z?   -  person    schedule 09.03.2014
comment
Каждый файл содержит 00z, 06z, 12z и 18z (время, UTC). это данные 4 раза в день. так, скажем, в файле t2m.2000.nc, t=1464, 4 раза в день в течение одного года. Поскольку данные находятся на поверхности (температура 2 метра), значения Z = 1. Это глобальные данные с привязкой к сетке.   -  person Fir Nor    schedule 10.03.2014


Ответы (1)


Итак, t2mtr - это трехмерный массив.

ntimes, ny, nx = shape(t2mtr)

Это суммирует все значения по 1-й оси:

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

Другим вариантом было бы изменить форму tm2tr (не работает для високосных лет).

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

Вы можете проверить выборку времени с помощью чего-то вроде:

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

Разве в наборе данных нет временной переменной? Возможно, вы сможете извлечь из этого нужные временные индексы. Я работал с данными ЕЦСПП несколько лет назад, но не помню многих деталей.

Что касается вашей ошибки contourf, я бы проверил форму трех основных аргументов: x, y, t2mtr. Они должны совпадать. Я не работал с Basemap.

person hpaulj    schedule 09.03.2014
comment
Спасибо @hpaulj за объяснение и решение. Именно то, что я ищу. Да, данные могут быть установлены в определенные индексы времени, но я скачал так, как сделал, чтобы сэкономить время (и немного места). Еще раз спасибо. - person Fir Nor; 10.03.2014
comment
привет @hpaulj, мне интересно в yr_starts = np.arange(0,3)*365 , относится ли (0,3) к временному шагу (в данном случае 4 часа в день) или к какому-то случайному числу лет. Спасибо! - person Fir Nor; 13.03.2014
comment
np.arange(0,3)*365 — это всего лишь [0, 365, 2*365], которые, как я предполагал, будут индексами последовательных точек данных за 1 января. Ему нужны дополнительные *6 для данных за 4 часа и смещение для других начальных точек и високосных лет. - person hpaulj; 13.03.2014
comment
Я бы также проверил данные для nan (или эквивалента). Когда я использовал ECMWF, мне нужны были как точечные данные, такие как температура, так и данные за период, такие как осадки. Один набор имел nan конец запрошенного периода времени, другой — его начало. - person hpaulj; 13.03.2014
comment
Спасибо @hpaulj .. Я только что понял, что также могу проверить консоль numpy.arangeon python. :) - person Fir Nor; 14.03.2014