Объединение большого количества файлов netCDF

У меня есть большая папка с файлами netCDF (.nc), каждый из которых имеет похожее имя. Файлы данных содержат переменные времени, долготы, широты и ежемесячных осадков. Цель состоит в том, чтобы получить среднемесячное количество осадков за X лет для каждого месяца. В итоге у меня будет 12 значений, представляющих среднемесячное количество осадков за X лет для каждой широты и долготы. Каждый файл находится в одном и том же месте на протяжении многих лет. Каждый файл начинается с того же имени и заканчивается «date.sub.nc», например:

'data1.somthing.somthing1.avg_2d_Ind_Nx.200109.SUB.nc'
'data1.somthing.somthing1.avg_2d_Ind_Nx.200509.SUB.nc'
'data2.somthing.somthing1.avg_2d_Ind_Nx.201104.SUB.nc'
'data2.somthing.somthing1.avg_2d_Ind_Nx.201004.SUB.nc'
'data2.somthing.somthing1.avg_2d_Ind_Nx.201003.SUB.nc'
'data2.somthing.somthing1.avg_2d_Ind_Nx.201103.SUB.nc'
'data1.somthing.somthing1.avg_2d_Ind_Nx.201203.SUB.nc'

Концовка - YearMonth.SUB.nc. Пока что у меня есть:

array=[]
f = nc.MFDataset('data*.nc')
precp = f.variables['prectot']
time = f.variables['time']
array = f.variables['time','longitude','latitude','prectot'] 

Я получаю KeyError: ('время', 'долгота', 'широта', 'пректот'). Есть ли способ объединить все эти данные, чтобы я мог ими управлять?


person BBHuggin    schedule 09.03.2015    source источник
comment
Что вы подразумеваете под объединением данных? Это уже все в одном объекте MFDataset благодаря вашей строке f = nc.MFDataset.... Другими словами, массив f.variables['prectot'][:] уже является трехмерным массивом с размерами (время, широта, долгота), содержащим значения prectot для каждого (время, широта, долгота).   -  person Spencer Hill    schedule 10.03.2015
comment
Кроме того, re: your KeyError, f.variables - это Dict, и для любого Dict вы можете получить доступ только к одному из его значений за раз, то есть f.variables['time'] или f.variables['longitude'], а не f.variables['time', 'longitude']. Но, как было сказано в моем предыдущем комментарии, все, что вам нужно, это f.variables['prectot'] (при условии, что я правильно вас понимаю).   -  person Spencer Hill    schedule 10.03.2015
comment
Понятно, я не был уверен в том, что на самом деле делает MFDataset. Я попробовал функцию glob.glob, но она просто составила список всех моих файлов. Спасибо.   -  person BBHuggin    schedule 10.03.2015


Ответы (3)


Как упомянул @CharlieZender, ncra - это то, что вам нужно, и я дам более подробную информацию об интеграции этой функции в скрипт Python. (PS - вы можете легко установить NCO с помощью Homebrew, например http://alejandrosoto.net/blog/2014/01/22/setting-up-my-mac-for-scientific-research/)

import subprocess
import netCDF4
import glob
import numpy as np

for month in range(1,13):
    # Gather all the files for this month
    month_files = glob.glob('/path/to/files/*{0:0>2d}.SUB.nc'.format(month))


    # Using NCO functions ---------------
    avg_file = './precip_avg_{0:0>2d}.nc'.format(month)

    # Concatenate the files using ncrcat
    subprocess.call(['ncrcat'] + month_files + ['-O', avg_file])

    # Take the time (record) average using ncra 
    subprocess.call(['ncra', avg_file, '-O', avg_file])

    # Read in the monthly precip climatology file and do whatever now
    ncfile = netCDF4.Dataset(avg_file, 'r')
    pr = ncfile.variables['prectot'][:,:,:]
    ....

    # Using only Python -------------
    # Initialize an array to store monthly-mean precip for all years
    # let's presume we know the lat and lon dimensions (nlat, nlon)
    nyears = len(month_files)
    pr_arr = np.zeros([nyears,nlat,nlon], dtype='f4')

    # Populate pr_arr with each file's monthly-mean precip
    for idx, filename in enumerate(month_files):
        ncfile = netCDF4.Dataset(filename, 'r')
        pr = ncfile.variable['prectot'][:,:,:]  
        pr_arr[idx,:,:] = np.mean(pr, axis=0)
        ncfile.close()

    # Take the average along all years for a monthly climatology
    pr_clim = np.mean(pr_arr, axis=0)  # 2D now [lat,lon]
person N1B4    schedule 10.03.2015
comment
Вы можете использовать чистый Python для достижения того же результата, NCO просто очень помогает уменьшить требования к времени / памяти. Я отредактирую свой ответ, чтобы показать вам, как использовать только Python. - person N1B4; 12.03.2015
comment
Две вещи: во-первых, цикл for для месяца в диапазоне (1,13): принимает файлы только за декабрь. При изменении диапазона с 1,13 на 1,12 используются только файлы за ноябрь. Может ли другой цикл for над первым решить эту проблему? Во-вторых, как только я попадаю в цикл заполнения pr_arr for, ему не нравится pr_arr [idx,:,:] = np.mean (pr, axis = 0), говорящее ValueError: не удалось передать массив ввода из shape (5,5) в форму (0,0). Другой вопрос: что на самом деле делает {0: 0 ›2d}? - person BBHuggin; 13.03.2015
comment
Странно, что цикл for не работает, он работает с моей стороны. Сначала он должен выбрать все файлы января (месяц = ​​1), а затем перейти к следующему месяцу. Попробуйте распечатать информацию после циклов for, чтобы помочь диагностировать проблему. Какая форма prectot для начала? Убедитесь, что вы вставляете размеры nlat и nlon в инициализацию pr_arr. '{0: 0 ›2d} дополняет целые числа начальным 0, например 1 становится 01, именно так структурированы имена файлов. Ознакомьтесь с дополнительной информацией о форматировании строк для получения дополнительных сведений об этом. - person N1B4; 13.03.2015
comment
Исправил проблему с петлей. Я добавил функцию добавления, теперь она выглядит так. month_files.append (glob.glob (* {0: 0 ›2d} .SUB.nc'.format (месяц))). Я видел, что при печати он печатал все файлы с 1-го по 12-й месяцы, но, похоже, вводил только последний набор (декабрь-12) в month_files. Когда вы говорите «Предположим, что широта и долгота» означает, что мне нужно создать новый массив с такими же широтой и долготой в моих файлах .nc? prectot имеет форму (396L, 5L, 5L). - person BBHuggin; 13.03.2015
comment
Вам просто нужно вставить 5, где nlat и nlon находятся в строке, где определено pr_arr. Кстати, если файлы хранятся ежемесячно - для начала, тогда форма prectot должна быть (1 x шир x дол) ... Я не уверен, откуда взялся 396. В этом случае вам не нужно вычислять среднемесячное значение, как я показал выше. Просто заполните pr_arr[idx,:,;] среднемесячным осадком из файла за каждый месяц, т.е. pr_arr[idx,:,:] = pr[0,:,:] - person N1B4; 14.03.2015
comment
396 - это данные за 33 года. (т.е. 396/12 = 33). Строка с [ncfile = netCDF4.Dataset (filename, 'r')] имеет проблему. Выдает «TypeError: ожидаемая строка или объект Unicode, список найден». - person BBHuggin; 15.03.2015

NCO делает это с

ncra *.01.SUB.nc pcp_avg_01.nc
ncra *.02.SUB.nc pcp_avg_02.nc
...
ncra *.12.SUB.nc pcp_avg_12.nc
ncrcat pcp_avg_??.nc pcp_avg.nc

Конечно, первые двенадцать команд могут быть выполнены с помощью цикла Bash, уменьшив общее количество строк до менее пяти. Если вы предпочитаете писать скрипты на Python, вы можете проверить свои ответы с помощью этого. ncra находится здесь.

person Charlie Zender    schedule 10.03.2015
comment
Я использую Enthought Canopy для python, и мне не удалось найти или загрузить пакет NCO для получения функции ncra. Я бы хотел его использовать. - person BBHuggin; 10.03.2015
comment
Думаю, я установил NCO. когда я запускаю ncra * 01.SUB.nc pcp_avg_01.nc, он выводит {SyntaxError: недопустимый синтаксис ncra * 01.SUB.nc pcp_avg_01.nc} с курсором, указывающим на S в SUB.nc, не знаю, как это исправить. - person BBHuggin; 16.03.2015
comment
Вы работаете в Windows, поэтому подстановка оболочки не работает. Замените * 01.SUB.nc фактическим списком файлов для ввода. Руководство может быть полезным. См. nco.sourceforge.net/nco.html#Specifying-Input-Files - person Charlie Zender; 17.03.2015

Команда ymonmean вычисляет среднее значение календарных месяцев в CDO. Таким образом, задачу можно выполнить двумя строками:

cdo mergetime data*.SUB.nc  merged.nc  # put files together into one series
cdo ymonmean merged.nc annual_cycle.nc # mean of all Jan,Feb etc. 

cdo может также рассчитать годовой цикл другой статистики, ymonstd, ymonmax и т. д., а единицы времени могут быть днями или пентадами, а также месяцами. (например, ydaymean).

person Adrian Tompkins    schedule 06.09.2017