Конвертируем netCDF в ASCII при помощи Python в Windows

Николай Колдунов
[email protected]

koldunov.net

Задача: Помочь друзьям виндузятникам сконвертировать netCDF в ASCII, попутно установив на их компьютеры Python, в надежде, что они таки постепенно забудут про дельфи, фортран и прочие гадости. Заодно попробовать удобно ли в ipython notebook писать посты.

Инструменты: cdo, Pyhton(x,y), ipython notebook

Люди работающие под виндоуз любят ASCII, я знаю, я сам был такой. Они готовы переводить в ASCII всё на свете, включая данные моделей IPCC, которые занимая в бинарном формате сотни гигабайт, будучи переведены в ASCII превращаются в монстров, сжирающих всё доступное дисковое пространство в радиусе нескольких километров. Но таковы реалии виндузовой жизни, многие программы там, особенно занимающиеся отрисовкой и анализом данных, хотят, чтобы им скармливали текст и только текст.

У проблемы перевода из netCDF в ASCII существует множество решений. Можно сделать дамп всего файла, заголовка, или отдельной переменной при помощи программки ncdump.exe. Небольшую инструкцию как это сделать и где взять эту неуловимую программку можно почитать здесь (там пишут про HDF, но для netCDF эта инструкция также подходит). Правда этот дамп вам придётся потом ещё долго и печально разбирать, поскольку то, что вы увидите, будет довольно сильно отличаться от желаемой таблички время/широта/долгота/значение (я знаю мечтаете вы именно об этом :)).

Тем кто знаком с Matlab с некоторых пор вообще стало хорошо, поскольку там появилась поддержка netCDF файлов из коробки и о том, как их там открывать можно почитать тут.

Здесь я расскажу как сконвертировать netCDF в ASCII при помощи Python, при этом формат вывода вы сможете задавать какой пожелаете. Упражняться будем, как обычно, на файлах NCEP реанализа.

Установка Pyhton(x,y) и cdo

Учёному, который хочет писать скрипты на питон под Windows, нужно кроме самого питона установить ещё много полезных в хозяйстве модулей, таких как numpy, scipy, ipyhton. Дело это зачастую муторное, особенно для такого новичка в виндоуз администрировании как я. Однако нашлись добрые люди, собрали всё необходимое в один большой пакет и распространяют его под именем Pyhton(x,y). Преимущество в том, что не нужно заморачиваться с самостоятельной настройкой пакетов, а недостаток, на мой взгляд, в том, что получившийся экзешник весит больше 500 мегабайт, что, к сожалению, пока ещё не тот объём, который можно безболезненно загрузить в каждом Нигерийском селе.

Процесс установки после скачивания прост и стандартен. На каком-то этапе вам будет предложено выбрать нужные пакеты - ставьте всё, что помечено. При том, что установится у вас куча всего, нам этого мало :) Чтобы рисовать красивые карты, нам необходим пакет Basemap, который устанавливается как дополнительный модуль к Pyhton(x,y), и лежит тут (первый в списке). Установив его вы получите на вашей Win машине полноценную среду для научного программирования под питон.

На самом деле я с трудом допускаю мысль, что вы хотите переводить в текстовый формат, скажем весь файл с реанализом ежедневных значений приземной температуры за год. Было бы неплохо сначала проделать над ним какие-то операции - осреднить по месяцам, или сезонам, построить стандартное отклонение наконец. Для этого человечество пока не придумало ничего лучше, чем climate data operators, или cdo. Небольшой рассказ о том, зачем они нужны и как ими пользоваться можно, как обычно найти на koldunov.net. С недавних пор они появились и под Win, так что скачивайте последнюю версию отсюда и распаковывайте её в ту папку, где вы будете работать с netCDF. В моём случае это папка C:/notebook/ (естественно вам её нужно будет предварительно создать).

Запуск IPython notebook

Ещё один маленький шажок и мы сможем перейти к работе над netCDF. Работать мы будем не из питоновской консоли, и не станем писать скрипты в отдельных файлах, а потом запускать их. Мы будем делать всё в IPython notebook, удивительном изобретении ребят из IPython, в котором я пишу этот пост, и который позволит вам одновременно запускать код и читать описание происходящего. Сперва загрузите ноутбук из которого сделан этот пост отсюда в вашу рабочую папку (C:/notebook/). Затем откройте консоль (если вы на XP то это Пуск >> Программы >> Стандартные >> Командная строка) и перейдите в свою рабочую папку (используйте команду cd). Если ваша папка, также как у меня, расположена по адресу C:/notebook/ то вам нужно будет ввести:

cd \
cd notebook

И наконец запускаем IPyhton notebook, выполнив находясь в вашем рабочем каталоге команду

ipython notebook

В вашем броузере должна открыться вкладка со списком ноутбуков, в котором, скорее всего, будет находиться только ноутбук, который вы загрузили, и, надеюсь, уже сейчас читаете. Теперь, когда все приготовления закончены, давайте перейдём к работе с нашими любимыми netCDF файлами.

Подготовка к экспорту

Для начала загружаем файл NCEP реанализа в рабочую папку. Я буду использовать файл ежедневной приземной температуры воздуха за 2012 год. IPython позволяет делать системные вызовы (перед ними должен стоять восклицательный знак), так что мы можем посмотреть, что у нас в папке:

In [3]:
!dir

 ’®¬ ў гбва®©б⢥ C ­Ґ Ё¬ҐҐв ¬ҐвЄЁ.
 ‘ҐаЁ©­л© ­®¬Ґа ⮬ : 9CDD-02A9

 ‘®¤Ґа¦Ё¬®Ґ Ї ЇЄЁ C:\notebook

15.02.2013  21:08    <DIR>          .
15.02.2013  21:08    <DIR>          ..
15.02.2013  21:08         7я700я344 air.sig995.2012.nc
07.01.2013  11:26        10я192я583 cdo.exe
15.02.2013  21:07            29я927 netCDF_to_ASCII_with_Python.ipynb
06.08.2012  13:40           127я720 pthreadGC2.dll
               4 д ©«®ў     18я050я574 Ў ©в
               2 Ї Ї®Є   2я759я348я224 Ў ©в бў®Ў®¤­®

Понятно, что у нас в папке находится экзешник cdo, дополнительная библиотека неоходимая ему для работы, файл с ноутбуком и, собственно наш netCDF файл. Посмотрим, что нам расскажут об этом файле cdo

In [11]:
!cdo pardes air.sig995.2012.nc

  -1  air           mean Daily Air temperature at sigma level 995 [degK]

cdo pardes: Processed 1 variable
In [12]:
!cdo showyear air.sig995.2012.nc

 2012

cdo showyear: Processed 1 variable over 366 timesteps
In [13]:
!cdo showmon air.sig995.2012.nc

  1  2  3  4  5  6  7  8  9 10 11 12

cdo showmon: Processed 1 variable over 366 timesteps

Мы узнали имя переменной, а также какие в файле содержатся годы и месяцы. Теперь давайте посчитаем среднюю температуру за сезон:

In [15]:
!cdo seasmean air.sig995.2012.nc air.sig995.2012_sm.nc

cdo seasmean: Processed 3847392 values from 1 variable over 366 timesteps
In [16]:
!cdo sinfo air.sig995.2012_sm.nc

   File format: netCDF
    -1 : Institut Source   Ttype    Levels Num  Gridsize Num Dtype : Parameter ID
     1 : unknown  unknown  instant       1   1     10512   1  I16  : -1         
   Grid coordinates :
     1 : lonlat       > size      : dim = 10512  nlon = 144  nlat = 73
                        lon       : first = 0  last = 357.5  inc = 2.5  degrees_east  circular
                        lat       : first = 90  last = -90  inc = -2.5  degrees_north
   Vertical coordinates :
     1 : surface                  : 0 
   Time coordinate :  5 steps
     RefTime =  0001-01-01 00:00:00  Units = hours  Calendar = STANDARD
  YYYY-MM-DD hh:mm:ss  YYYY-MM-DD hh:mm:ss  YYYY-MM-DD hh:mm:ss  YYYY-MM-DD hh:mm:ss
  2012-02-29 00:00:00  2012-05-31 00:00:00  2012-08-31 00:00:00  2012-11-30 00:00:00
  2012-12-31 00:00:00

cdo sinfo: Processed 1 variable over 5 timesteps

Как видите получилось у нас не четыре, а пять полей, поскольку зим за 2012 год было аж две штуки.

Допустим, весь мир нам не нужен, а мы хотим вырезать только Европу:

In [20]:
!cdo sellonlatbox,0,90,30,70 air.sig995.2012_sm.nc air.sig995.2012_sm_europe.nc

cdo sellonlatbox: Processed 52560 values from 1 variable over 5 timesteps

Этот файл и будем конвертировать в ASCII.

Открываем netCDF файл в Pyhton

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

In [19]:
from scipy.io import netcdf
import numpy

Теперь открываем файл:

In [21]:
f = netcdf.netcdf_file('air.sig995.2012_sm_europe.nc')

Смотрим, что там у нас внутри:

In [22]:
f.variables
Out[22]:
{'air': <scipy.io.netcdf.netcdf_variable at 0x2c43bb0>,
 'lat': <scipy.io.netcdf.netcdf_variable at 0x2c436b0>,
 'lon': <scipy.io.netcdf.netcdf_variable at 0x2c43690>,
 'time': <scipy.io.netcdf.netcdf_variable at 0x2c43870>}

Нам понадобятся все эти переменные, так что импортируем их в переменны питона:

In [23]:
air = f.variables['air']
lat = f.variables['lat']
lon = f.variables['lon']
time = f.variables['time']

Если вы теперь поставите точку после имени переменной (любой, не только air) и нажмёте TAB, то увидите список различных атрибутов этой переменной, предлагаю вам их поизучать:

In [32]:
air.

Два из них особенно важны, это add_offset и scale_factor. Откуда они взялись и зачем нужны можно почитать здесь. Если читать лень, то просто знайте: чтобы получить нормальные данные в Кельвинах нужно умножить air на scale_factor и прибавить add_offset. Так надо.

В данный момент в наших переменных содержатся не только данные, но и дополнительная информация, вроде атрибутов. Для того, чтобы работать только с данным, мы можем их скопировать в отдельную переменную, попутно сконвертировав их в Кельвины:

In [35]:
airK = (air.data*air.scale_factor)+air.add_offset

Ещё вариант

In [39]:
airK = (air[:]*air.scale_factor)+air.add_offset

В принципе, если бы мы не хотели получать информацию об атрибутах переменных, то могли бы импортировать только данные следующим образом:

air = f.variables['air'][:]
lat = f.variables['lat'][:]
lon = f.variables['lon'][:]
time = f.variables['time'][:]

но в этом случае информацию, например, о scale_factor и add_offset пришлось бы добывать из других источников.

Переведём данные из Кельвинов в Цельсии

In [41]:
airC = airK-273.15

Вы уже, наверное, поняли, что никаких циклов с перебиранием каждого значения и вычитанием магической цифры тут не нужно, всё как в матлабе, операция вычитания производится для всех элементов массива.

Теперь взгляним на наши координаты:

In [42]:
lat.shape
Out[42]:
(17,)
In [43]:
lon.shape
Out[43]:
(37,)

Они одномерные. Ничего плохого в этом нет, но мы хотим нарисовать наши данные, прежде чем экспортировать, а для этого, нужно перевсти широты и долготы в двумерный формат (то есть двумерному массиву, в котором содержится поле с данными, будут соответствовать два двумерных массива в широтами и долготами).

In [45]:
lons, lats = numpy.meshgrid(lon[:],lat[:])

Поглядим, в каком формате у нас время

In [70]:
time.units
Out[70]:
'hours since 1-01-01 00:00:00'

Время у нас в часах от начала нашей эры. Выглядят они вот так:

In [71]:
time[:][0]
Out[71]:
17629512.0

Чтобы перевести эти часы в удобоваримый формат, нам нужно воспользоваться модулем datetime

In [72]:
import datetime

Создадим объект специального типа datetime, при помощи которого удобно управляться со временем, в котором будет дата с которой считаются часы в нашем формате времени:

In [154]:
first_date = datetime.datetime(1,1,1,0,0,0)
In [155]:
type(first_date)
Out[155]:
datetime.datetime

Теперь представим нашу первую дату в виде временной разницы

In [168]:
difference1 = datetime.timedelta(hours=(time[:][0]-48))

По причине, в которой я ещё не разобрался, эта функция ошибается на два дня, так что мы их (то есть 48 часов) вычитаем.

Прибавим к начальной дате разницу, и получим datetime переменную, в которой будет содержаться наша дата.

In [169]:
present = first_date+difference1

У переменной present много интересных методов, которые позволяют узнать не только дату, но и, например, день недели. Но нас будет интересовать полная запись даты

In [170]:
present.ctime()
Out[170]:
'Wed Feb 29 00:00:00 2012'

Как видите в качестве даты даётся последний день сезона. Информацию о времени мы используем позже при отрисовке данных и экспорте в netCDF.

Теперь немножно магии, чтобы отображать картинки внутри ноутбука:

In [46]:
%pylab inline

Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.

Импортируем модуль для отрисовки карт:

In [47]:
from mpl_toolkits.basemap import Basemap

Генерируем карту

In [65]:
m = Basemap(llcrnrlon=-10,llcrnrlat=25,urcrnrlon=95,urcrnrlat=75,projection='mill', resolution='l')

Переводим широты и долготы в координаты карты

In [66]:
x, y = m(lons, lats)

И отрисовываем наши данные:

In [171]:
# Создаём рисунок и задаём его размеры
fig = plt.figure(figsize=(10,10))
# Выбираем сезон
season=2
# Рисуем береговую линию
m.drawcoastlines(linewidth=1)
# Рисуем границы стран
m.drawcountries(linewidth=0.5)
# Рисуем границу карты
m.drawmapboundary
# Создаём карту с залитыми изолинииями, 20 интервалов
cs = m.contourf(x,y,airC[season,:,:],20)
# Создаём легенду
cbar = m.colorbar(cs,location='bottom',pad="5%")
# Подписываем легенду
cbar.set_label(r'$^{\circ}\mathrm{C}$')
# Переводим время в человеческое
difference1 = datetime.timedelta(hours=(time[:][season]-48))
present = first_date+difference1
# Отрисовываем заголовок
plt.title(present.ctime())
Out[171]:
<matplotlib.text.Text at 0x3b31b70>

С определением Европы мы явно промахнулись :) Чтобы исправить эту ситуацию можно вернуться к шагу, в котором мы вырезаем регион и подправить координаты, которые передаются в cdo.

В принципе, загрузив данные в питон (назовём это так) вы можете делать с ними всё, что угодно. В этом смысле процесс не сильно отличается от работы, например, в Matlab. Но мы с вами собирались сохранить данные в ASCII, давайте этим и займёмся.

Экспорт в ASCII

Алгоритм сохранения прост: открываем файл, последовательно записываем в него строчки в нужном нам формате, закрываем файл. Теперь посмотрим как это выглядит на практике. Сначала пример, в котором каждое поле сохраняется в отдельный файл:

In [221]:
#Цикл по времени
for one_field in range(airC.shape[0]):
    
    #Открываем файл, дописываем к имени номер поля и нули, 
    #чтобы при большом количестве файлов они нормально сортировались
    ofile = open("netcdf_to_ASCII_field_"+str(one_field).zfill(5)+".txt","w")
    
    #Переводим время в человесеский формат
    difference1 = datetime.timedelta(hours=(time[:][one_field]-48))
    present = first_date+difference1
    
    #Печатаем, чтобы видеть прогресс
    print(present.ctime())
    
    #Цикл по строкам
    for nn in range(airC.shape[1]):
        
        #Цикл по столбцам
        for kk in range(airC.shape[2]):
            
            #Запись строки в файл. Сначала перечисляем форматы,
            #в которых будт выводиться переменные, а потом перечисляем
            #сами переменные
            ofile.write('%5.2f  %5.2f %5.2f %5.0f %5.0f %5.0f\n' % \
            (lons[nn,kk], lats[nn,kk], airC[one_field,nn,kk], present.day, present.month, present.year))
                        
                        
    ofile.close()
            

Wed Feb 29 00:00:00 2012
Thu May 31 00:00:00 2012
Fri Aug 31 00:00:00 2012
Fri Nov 30 00:00:00 2012
Mon Dec 31 00:00:00 2012

Самое "сложное", или скорее неочевидное тут это форматы вывода данных. Про них можно почитать тут.

Вариант с выводом в один файл с дата стампами между полями:

In [224]:
#Открывать будем только один файл
ofile = open("netcdf_to_ASCII_all_fields.txt","w")

#Цикл по времени
for one_field in range(airC.shape[0]):
        
    #Переводим время в человесеский формат
    difference1 = datetime.timedelta(hours=(time[:][one_field]-48))
    present = first_date+difference1
    
    #Печатаем, чтобы видеть прогресс
    print(present.ctime())
    
    #Записываем время в начале поля,
    #тут мы решили ещё и часы с минутами приписать
    ofile.write('%5.0f  %5.0f %5.0f %5.0f %5.0f\n' % (present.day, present.month, present.year, present.hour, present.minute ))
    
    #Цикл по строкам
    for nn in range(airC.shape[1]):
        
        #Цикл по столбцам
        for kk in range(airC.shape[2]):
            
            #Запись строки в файл. Сначала перечисляем форматы,
            #в которых будт выводиться переменные, а потом перечисляем
            #сами переменные
            ofile.write('%5.2f  %5.2f %5.2f\n' % \
            (lons[nn,kk], lats[nn,kk], airC[one_field,nn,kk]))
                        
                        
ofile.close()

Wed Feb 29 00:00:00 2012
Thu May 31 00:00:00 2012
Fri Aug 31 00:00:00 2012
Fri Nov 30 00:00:00 2012
Mon Dec 31 00:00:00 2012

Как видите слагаемые буквально чуть-чуть поменялись местами, а результат изменился :)

Пользуясь двумя этими шаблонами, как мне кажется, можно сохранить большинство данных, которые вам могут понадобиться, в текстовых форматах, которые нужны именно вам, а не каким-то совершенно чужим людям :) Однако научившись "загружать" и отображать данные в питоне, я надеюсь вам захочется провести в нём больше времени, и может в один прекрасный день вы поймёте, что ASCII файлы вам больше не нужны :)