Николай Колдунов
koldunovn@gmail.com
Задача: Помочь друзьям виндузятникам сконвертировать 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 реанализа.
Учёному, который хочет писать скрипты на питон под Windows, нужно кроме самого питона установить ещё много полезных в хозяйстве модулей, таких как numpy, scipy, ipyhton. Дело это зачастую муторное, особенно для такого новичка в виндоуз администрировании как я. Однако нашлись добрые люди, собрали всё необходимое в один большой пакет и распространяют его под именем Pyhton(x,y). Преимущество в том, что не нужно заморачиваться с самостоятельной настройкой пакетов, а недостаток, на мой взгляд, в том, что получившийся экзешник весит больше 500 мегабайт, что, к сожалению, пока ещё не тот объём, который можно безболезненно загрузить в каждом Нигерийском селе.
Процесс установки после скачивания прост и стандартен. На каком-то этапе вам будет предложено выбрать нужные пакеты - ставьте всё, что помечено. При том, что установится у вас куча всего, нам этого мало :) Чтобы рисовать красивые карты, нам необходим пакет Basemap, который устанавливается как дополнительный модуль к Pyhton(x,y), и лежит тут (первый в списке). Установив его вы получите на вашей Win машине полноценную среду для научного программирования под питон.
На самом деле я с трудом допускаю мысль, что вы хотите переводить в текстовый формат, скажем весь файл с реанализом ежедневных значений приземной температуры за год. Было бы неплохо сначала проделать над ним какие-то операции - осреднить по месяцам, или сезонам, построить стандартное отклонение наконец. Для этого человечество пока не придумало ничего лучше, чем climate data operators, или cdo. Небольшой рассказ о том, зачем они нужны и как ими пользоваться можно, как обычно найти на koldunov.net. С недавних пор они появились и под Win, так что скачивайте последнюю версию отсюда и распаковывайте её в ту папку, где вы будете работать с netCDF. В моём случае это папка C:/notebook/ (естественно вам её нужно будет предварительно создать).
Ещё один маленький шажок и мы сможем перейти к работе над netCDF. Работать мы будем не из питоновской консоли, и не станем писать скрипты в отдельных файлах, а потом запускать их. Мы будем делать всё в IPython notebook, удивительном изобретении ребят из IPython, в котором я пишу этот пост, и который позволит вам одновременно запускать код и читать описание происходящего. Сперва загрузите ноутбук из которого сделан этот пост отсюда в вашу рабочую папку (C:/notebook/). Затем откройте консоль (если вы на XP то это Пуск >> Программы >> Стандартные >> Командная строка) и перейдите в свою рабочую папку (используйте команду cd). Если ваша папка, также как у меня, расположена по адресу C:/notebook/ то вам нужно будет ввести:
cd \
cd notebook
И наконец запускаем IPyhton notebook, выполнив находясь в вашем рабочем каталоге команду
ipython notebook
В вашем броузере должна открыться вкладка со списком ноутбуков, в котором, скорее всего, будет находиться только ноутбук, который вы загрузили, и, надеюсь, уже сейчас читаете. Теперь, когда все приготовления закончены, давайте перейдём к работе с нашими любимыми netCDF файлами.
Для начала загружаем файл NCEP реанализа в рабочую папку. Я буду использовать файл ежедневной приземной температуры воздуха за 2012 год. IPython позволяет делать системные вызовы (перед ними должен стоять восклицательный знак), так что мы можем посмотреть, что у нас в папке:
!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
!cdo pardes air.sig995.2012.nc
-1 air mean Daily Air temperature at sigma level 995 [degK]
cdo pardes: Processed 1 variable
!cdo showyear air.sig995.2012.nc
2012
cdo showyear: Processed 1 variable over 366 timesteps
!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
Мы узнали имя переменной, а также какие в файле содержатся годы и месяцы. Теперь давайте посчитаем среднюю температуру за сезон:
!cdo seasmean air.sig995.2012.nc air.sig995.2012_sm.nc
cdo seasmean: Processed 3847392 values from 1 variable over 366 timesteps
!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 год было аж две штуки.
Допустим, весь мир нам не нужен, а мы хотим вырезать только Европу:
!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, нам нужно импортировать пару модулей:
from scipy.io import netcdf
import numpy
Теперь открываем файл:
f = netcdf.netcdf_file('air.sig995.2012_sm_europe.nc')
Смотрим, что там у нас внутри:
f.variables
{'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>}
Нам понадобятся все эти переменные, так что импортируем их в переменны питона:
air = f.variables['air']
lat = f.variables['lat']
lon = f.variables['lon']
time = f.variables['time']
Если вы теперь поставите точку после имени переменной (любой, не только air) и нажмёте TAB, то увидите список различных атрибутов этой переменной, предлагаю вам их поизучать:
air.
Два из них особенно важны, это add_offset и scale_factor. Откуда они взялись и зачем нужны можно почитать здесь. Если читать лень, то просто знайте: чтобы получить нормальные данные в Кельвинах нужно умножить air на scale_factor и прибавить add_offset. Так надо.
В данный момент в наших переменных содержатся не только данные, но и дополнительная информация, вроде атрибутов. Для того, чтобы работать только с данным, мы можем их скопировать в отдельную переменную, попутно сконвертировав их в Кельвины:
airK = (air.data*air.scale_factor)+air.add_offset
Ещё вариант
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 пришлось бы добывать из других источников.
Переведём данные из Кельвинов в Цельсии
airC = airK-273.15
Вы уже, наверное, поняли, что никаких циклов с перебиранием каждого значения и вычитанием магической цифры тут не нужно, всё как в матлабе, операция вычитания производится для всех элементов массива.
Теперь взгляним на наши координаты:
lat.shape
(17,)
lon.shape
(37,)
Они одномерные. Ничего плохого в этом нет, но мы хотим нарисовать наши данные, прежде чем экспортировать, а для этого, нужно перевсти широты и долготы в двумерный формат (то есть двумерному массиву, в котором содержится поле с данными, будут соответствовать два двумерных массива в широтами и долготами).
lons, lats = numpy.meshgrid(lon[:],lat[:])
Поглядим, в каком формате у нас время
time.units
'hours since 1-01-01 00:00:00'
Время у нас в часах от начала нашей эры. Выглядят они вот так:
time[:][0]
17629512.0
Чтобы перевести эти часы в удобоваримый формат, нам нужно воспользоваться модулем datetime
import datetime
Создадим объект специального типа datetime, при помощи которого удобно управляться со временем, в котором будет дата с которой считаются часы в нашем формате времени:
first_date = datetime.datetime(1,1,1,0,0,0)
type(first_date)
datetime.datetime
Теперь представим нашу первую дату в виде временной разницы
difference1 = datetime.timedelta(hours=(time[:][0]-48))
По причине, в которой я ещё не разобрался, эта функция ошибается на два дня, так что мы их (то есть 48 часов) вычитаем.
Прибавим к начальной дате разницу, и получим datetime переменную, в которой будет содержаться наша дата.
present = first_date+difference1
У переменной present много интересных методов, которые позволяют узнать не только дату, но и, например, день недели. Но нас будет интересовать полная запись даты
present.ctime()
'Wed Feb 29 00:00:00 2012'
Как видите в качестве даты даётся последний день сезона. Информацию о времени мы используем позже при отрисовке данных и экспорте в netCDF.
Теперь немножно магии, чтобы отображать картинки внутри ноутбука:
%pylab inline
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.
Импортируем модуль для отрисовки карт:
from mpl_toolkits.basemap import Basemap
Генерируем карту
m = Basemap(llcrnrlon=-10,llcrnrlat=25,urcrnrlon=95,urcrnrlat=75,projection='mill', resolution='l')
Переводим широты и долготы в координаты карты
x, y = m(lons, lats)
И отрисовываем наши данные:
# Создаём рисунок и задаём его размеры
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())
<matplotlib.text.Text at 0x3b31b70>
С определением Европы мы явно промахнулись :) Чтобы исправить эту ситуацию можно вернуться к шагу, в котором мы вырезаем регион и подправить координаты, которые передаются в cdo.
В принципе, загрузив данные в питон (назовём это так) вы можете делать с ними всё, что угодно. В этом смысле процесс не сильно отличается от работы, например, в Matlab. Но мы с вами собирались сохранить данные в ASCII, давайте этим и займёмся.
Алгоритм сохранения прост: открываем файл, последовательно записываем в него строчки в нужном нам формате, закрываем файл. Теперь посмотрим как это выглядит на практике. Сначала пример, в котором каждое поле сохраняется в отдельный файл:
#Цикл по времени
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
Самое "сложное", или скорее неочевидное тут это форматы вывода данных. Про них можно почитать тут.
Вариант с выводом в один файл с дата стампами между полями:
#Открывать будем только один файл
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 файлы вам больше не нужны :)