Конвертируем netCDF в ASCII при помощи Python в Windows
Задача: Помочь друзьям виндузятникам сконвертировать netCDF в ASCII, попутно установив на их компьютеры Python, в надежде, что они таки постепенно забудут про дельфи, фортран и прочие гадости. Заодно попробовать удобно ли в ipython notebook писать посты.
?нструменты: cdo, Pyhton(x,y), ipython notebook
Важно: Этот пост изначально целиком написан в ipython notebook, что значит вы сможете при наличии установленного ipython notebook (как это сделать под Windows смотри ниже) исполнять и изменять код представленный в посте. Читать пост с хорошим форматированием, так как оно выглядит в ноутбуке, а также загрузить файл ноутбука можно здесь (загрузка файла в правом верхнем углу). Также доступна pdf версия этого поста автоматически сгенерированная nbconverter. А ещё версия этого поста в виде мимимишной презентации reveal. А теперь собственно пост:
Люди работающие под виндоуз любят 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 notebook
? наконец запускаем IPyhton notebook, выполнив находясь в вашем рабочем каталоге команду
В вашем броузере должна открыться вкладка со списком ноутбуков, в котором, скорее всего, будет находиться только ноутбук, который вы загрузили, и, надеюсь, уже сейчас читаете. Теперь, когда все приготовления закончены, давайте перейдём к работе с нашими любимыми netCDF файлами.
Подготовка к экспорту
Для начала загружаем файл NCEP реанализа в рабочую папку. Я буду использовать файл ежедневной приземной температуры воздуха за 2012 год. IPython позволяет делать системные вызовы (перед ними должен стоять восклицательный знак), так что мы можем посмотреть, что у нас в папке:
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
Понятно, что у нас в папке находится экзешник cdo, дополнительная библиотека необходимая ему для работы, файл с ноутбуком и, собственно наш netCDF файл. Посмотрим, что нам расскажут об этом файле cdo
cdo pardes: Processed 1 variable
cdo showyear: Processed 1 variable over 366 timesteps
cdo showmon: Processed 1 variable over 366 timesteps
Мы узнали имя переменной, а также какие в файле содержатся годы и месяцы. Теперь давайте посчитаем среднюю температуру за сезон:
Мы узнали имя переменной, а также какие в файле содержатся годы и месяцы. Теперь давайте посчитаем среднюю температуру за сезон:
-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 год было аж две штуки.
Допустим, весь мир нам не нужен, а мы хотим вырезать только Европу:
Этот файл и будем конвертировать в ASCII.
Открываем netCDF файл в Pyhton
Чтобы иметь возможность открывать netCDF, нам нужно импортировать пару модулей:
import numpy
Теперь открываем файл:
Смотрим, что там у нас внутри:
'lat': <scipy.io.netcdf.netcdf_variable at 0x2c436b0>,
'lon': <scipy.io.netcdf.netcdf_variable at 0x2c43690>,
'time': <scipy.io.netcdf.netcdf_variable at 0x2c43870>}
Нам понадобятся все эти переменные, так что импортируем их в переменны питона:
lat = f.variables['lat']
lon = f.variables['lon']
time = f.variables['time']
Если вы теперь поставите точку после имени переменной (любой, не только air) и нажмёте TAB, то увидите список различных атрибутов этой переменной, предлагаю вам их поизучать:
Два из них особенно важны, это add_offset и scale_factor. Откуда они взялись и зачем нужны можно почитать здесь. Если читать лень, то просто знайте: чтобы получить нормальные данные в Кельвинах нужно умножить air на scale_factor и прибавить add_offset. Так надо.
В данный момент в наших переменных содержатся не только данные, но и дополнительная информация, вроде атрибутов. Для того, чтобы работать только с данным, мы можем их скопировать в отдельную переменную, попутно сконвертировав их в Кельвины:
Ещё вариант
В принципе, если бы мы не хотели получать информацию об атрибутах переменных, то могли бы импортировать только данные следующим образом:
lat = f.variables['lat'][:]
lon = f.variables['lon'][:]
time = f.variables['time'][:]
но в этом случае информацию, например, о scale_factor и add_offset пришлось бы добывать из других источников.
Переведём данные из Кельвинов в Цельсии
Вы уже, наверное, поняли, что никаких циклов с перебиранием каждого значения и вычитанием магической цифры тут не нужно, всё как в матлабе, операция вычитания производится для всех элементов массива.
Теперь взгляним на наши координаты:
Они одномерные. Ничего плохого в этом нет, но мы хотим нарисовать наши данные, прежде чем экспортировать, а для этого, нужно перевсти широты и долготы в двумерный формат (то есть двумерному массиву, в котором содержится поле с данными, будут соответствовать два двумерных массива в широтами и долготами).
Поглядим, в каком формате у нас время
Время у нас в часах от начала нашей эры. Выглядят они вот так:
Чтобы перевести эти часы в удобоваримый формат, нам нужно воспользоваться модулем datetime
Создадим объект специального типа datetime, при помощи которого удобно управляться со временем, в котором будет дата с которой считаются часы в нашем формате времени:
Теперь представим нашу первую дату в виде временной разницы
По причине, в которой я ещё не разобрался, эта функция ошибается на два дня, так что мы их (то есть 48 часов) вычитаем.
Прибавим к начальной дате разницу, и получим datetime переменную, в которой будет содержаться наша дата.
У переменной present много интересных методов, которые позволяют узнать не только дату, но и, например, день недели. Но нас будет интересовать полная запись даты
Как видите в качестве даты даётся последний день сезона. ?нформацию о времени мы используем позже при отрисовке данных и экспорте в netCDF.
Теперь немножно магии, чтобы отображать картинки внутри ноутбука:
For more information, type 'help(pylab)'.
?мпортируем модуль для отрисовки карт:
Генерируем карту
Переводим широты и долготы в координаты карты
? отрисовываем наши данные:
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())
С определением Европы мы явно промахнулись 🙂 Чтобы исправить эту ситуацию можно вернуться к шагу, в котором мы вырезаем регион и подправить координаты, которые передаются в cdo.
В принципе, загрузив данные в питон (назовём это так) вы можете делать с ними всё, что угодно. В этом смысле процесс не сильно отличается от работы, например, в Matlab. Но мы с вами собирались сохранить данные в ASCII, давайте этим и займёмся.
Экспорт в 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()
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()
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 файлы вам больше не нужны 🙂
Одно не понятно, зачем отдельно качать ipython, когда он уже есть в полном пакете python(x,y)?
@Oleg
А где написано, что его нужно отдельно качать?
Благодаря этому посту попробовал Python(x,y). Всё поставилось, но, iPython запускаться отказывался. Дело оказалось в непонятно откуда взявшихся mime-типах в реестре, в названиях которых была кириллица. Поправил это — и всё запустилось. Свой user experience описал здесь — http://thebaga.blogspot.com/2013/04/ipython.html Может, кто ещё с таким столкнётся.
Не знаю, где правильнее задать вопрос, но…А как можно прочитать файл, в котором есть пропущенные значения, на месте которых стоят просто пробелы? При этом в качестве delimeter выступает тоже пробел. ? еще одна загвоздка состоит в том, что иногда бывает несколько пропусков подряд, и тогда это может восприниматься за один..
Очень нужна помощь! В python еще совсем новичок…
2natavrame
Не могу сказать точно как вам помочь не видя самого файла ) Единственно, если количество символов в строке всегда одинаково, то можно открыв файл и считав из него строки разобрать каждую строку посимвольно. Например:
lines = ifile.readlines()
for line in lines:
param1 = line[0:5]
param2 = line[5:10]
Спасибо большое! Да, если задавать каждый параметр отдельно, то вроде нормально считывает. А дальше тогда надо поставить if-условие, что-то типа
if param==’\t\t\t\t\t\t’:
param=9999
Сам файл выглядит примерно так: (не знаю, как его можно загрузить сюда )
20674 1967 22.0 35.0 28.0 39.0 22.0 23.0 59.0 46.0 49.0 45.0 36.0 29.0
20674 1968 17.0 21.0 28.0 17.0 14.0 25.0 47.0 68.0 22.0 22.0
20674 1969 20.0 16.0 15.0 8.0 9.0 17.0 60.0 33.0 27.0 11.0 27.0 30.0
20674 1970 11.0 17.0 24.0 10.0 19.0 14.0 56.0 10.0 24.0 25.0 36.0 43.0
20674 1971 91.0 39.0 17.0 4.0 12.0 13.0 31.0 67.0 78.0 16.0 31.0 26.0
20674 1972 32.0 47.0 11.0 28.0 19.0 23.0 51.0 43.0 34.0 35.0 4.0 69.0
20674 1973 25.0 36.0 25.0 9.0 38.0 50.0 21.0 75.0 43.0 36.0 16.0 15.0
20674 1974 26.0 40.0 30.0 21.0 44.0 31.0 44.0 22.0 61.0 32.0 27.0
20674 1975 49.0 16.0 15.0 29.0 12.0 28.0 35.0 25.0 35.0 43.0 24.0 54.0
20674 1976 24.0 34.0 17.0 36.0 19.0 53.0 40.0 63.0 47.0 19.0 21.0 14.0
20674 1977
20674 1978 20.0 20.0 16.0 6.0 14.0 18.0 48.0 41.0 12.0 27.0 25.0 16.0
20674 1979 22.0 10.0 9.0 8.0 48.0 38.0 34.0 24.0 13.0 19.0 26.0
20674 1980 21.0 35.0 9.0 14.0 23.0 41.0 66.0 67.0 55.0 44.0 8.0 32.0
20674 1981 105.0 14.0 14.0 23.0 40.0 14.0 47.0 27.0 26.0 28.0
20674 1982 8.0 13.0 13.0 24.0 12.0 58.0 13.0 71.0 40.0 23.0 5.0 26.0
Это данные станционных наблюдений ВН??ГМ? по среднемесячным осадкам.
Да, по идее должно работать ) Если что обращайтесь )
Благодарю за помощь!!
Почему-то искать ’\t\t\t\t\t\t’ ему не понравилось и он решил оставлять все как есть. Поэтому мне помогли сделать задачу через фотран. Но, надеюсь, через некоторое время и у меня получится сделать это все-таки на python. )
2natavrame
Так мне стало грустно от того, что вам пришлось воспользоваться фортраном, что пришлось сделать ноутбук с объяснением как можно открыть ваши данные и начать их анализировать 🙂 Посмотреть его можно тут: http://nbviewer.ipython.org/5391791
Огромное Вам спасибо!!! Это просто волшебно! Чем больше читаю примеров использования, тем больше и больше нравится python. Особенно приятно удивило «Time series analysis with pandas», потому что то, что описывается в несколько строк в других языках, можно сделать всего лишь одной строчкой. ))
А можно задать еще один вопрос…может быть Вам приходилось сталкиваться с этим, можно посчитать площадь внутри замкнутых изолиний? ?ли это нудно описывать с самого нуля вручную?
*нужно
2natavrame
Я возможно не очень понял, что вы имеете в виду, но вот первое, что приходит в голову:
— если вы работаете с гридированными данными, то можно просто выбрать все гридбоксы с данными больше, меньше, или находящиеся между какими то значениями, и затем посчитать общую площадь, сложив площади всех удовлетворяющих условиям точек.
— если вы работаете с негридированными данными, то для того, чтобы получить изолинии вам их нужно будет интерполировать на какой-то грид, после чего см. пункт первый 🙂
Сделал маленький ноутбук с возможной реализацией этого дела http://nbviewer.ipython.org/5404107
нужно для дипломной работы…
подскажите пожалуйста, у меня данные координатной сетки не в долготе и широте, а в xx — абцисса и yy — ордината:
-25 xx abscissas
-26 yy ordinates
xx.data
Out[51]:
array([ 60000., 60000., 60000., 60000., 60000., 60000.,
60000., 60000., 60000., 60000., 60000., 60000.,
120000., 120000., 120000., 120000., 120000., 120000.,
120000., 120000., 120000., 120000., 120000., 120000.,
180000., 180000., 180000., 180000., 180000., 180000.,
180000., 180000., 180000., 180000., 180000., 180000.,
240000., 240000., 240000., 240000., 240000., 240000.,
240000., 240000., 240000., 240000., 240000., 240000.,
300000., 300000., 300000., 300000., 300000., 300000.,
300000., 300000., 300000., 300000., 300000., 300000.,
360000., 360000., 360000., 360000., 360000., 360000.,
360000., 360000., 360000., 360000., 360000., 360000.,
420000., 420000., 420000., 420000., 420000., 420000.,
420000., 420000., 420000., 420000., 420000., 420000.,
480000., 480000., 480000., 480000., 480000., 480000.,
480000., 480000., 480000., 480000., 480000., 480000.])
yy.data — выглядят аналогично
lon_0=0,lat_0=43
Как мне перевести абциссы с ординатами в широты с долготами или может строить выходя от имеющегося — то есть, может есть вариант как нанести данные на карту с такими вот хх и yy?
Подскажите, пожалуйста, как решить обратную задачу, т.е. из текстового файла соорудить netcdf на Рython.