Открываем netCDF в MATLAB (на примере NCEP реанализа)
Задача: Открыть файл формата netCDF в MATLAB
Решение: чистый незамутненный MATLAB
Постараюсь описать здесь ответ на этот животрепещущий вопрос, чтобы было куда отправлять страждущих с различных форумов. Начиная с версии 7.7, Matlab поддерживает работу с форматом netCDF нативно, без различных примочек, которые требовались раньше. Синтаксис, который используется для работы с netCDF файлами довольно странный, ну да, не мне судить. Здесь я опишу, как в Matlab совершить одно простое действие, а именно прочитать данные из файла.
Если вы собираетесь работать с netCDF, то для начала неплохо бы ознакомиться с тем, что из себя представляет этот формат. В этом случае дальнейшее не будет вам казаться абракадаброй. Но в общем случае можно обойтись и без этих уникальных знаний, если все, что вам нужно, это выковырять данные из файла и забыть о netCDF как о страшном сне. ?нформация о функциях Матлаба для работы с netCDF расположена тут.
Для начала, давайте скачаем подопытный файл. Это будет файл с шестичасовой приземной температурой любимого всеми гидрометеорологами NCEP реанализа. Я буду использовать вот этот файл за 2010 год (20 мегабайт). Чтобы открыть netCDF файл, нужно иметь представление о его структуре, о тех переменных которые в нем находятся и об аттрибутах этих переменных. Для этого служит функция ncdisp.
В результате вы получите табличку с данными по вашему файлу. По большому счету вас должно интересовать только то что находится в разделе Variables, то есть информация о переменных, содержащихся в вашем файле. В данном примере это переменные lat (широты), lon (долготы), time(время прошедшее от 1-1-1 00:00:0.0 в часах), air (температура воздуха в Кельвинах).
UPD: Говорят что не во всех релизах Матлаба есть функция ncdisp. Если у вас тот самый запущенный случай, то для виндоуз вы скачиваете бинарники программы ncdump отсюда, для Убунту Линукс и ему подобных пишите в терминале:
Далее в терминале виндоуз или линукс вводите команду
по идее, если вы работаете в Виндоуз, то программа ncdump.exe должна быть в папке с файлом. Эта команда выведет вам примерно ту же информацию что и ncdisp.
Сначала давайте попробуем добыть главное — температуру воздуха. Перво-наперво нужно открыть файл, для того чтобы Matlab мог с ним работать:
По каким-то своим внутренним причинам функции, которые запрашивают данные из переменных netCDF, работают не с их именами, а с их порядковыми номерами, поэтому нам нужно узнать каков порядковый номер переменной air:
varid =
3
Получился порядковый номер 3. Мы в общем и сами могли это посчитать, зная, что счет переменных в данном случае идет с нуля, а не с единицы.
Теперь, зная порядковый номер netCDF переменной, мы можем скопировать значения из нее в переменную Матлаба:
В итоге вы получите трехмерную матрицу [долгота, широта, время]. Казалось бы задача решена, но не стоит торопиться. Если вы взгляните на значения в переменной data, то они будут очень отдаленно напоминать то, что вы ожидаете увидеть, таких температур, даже если они выражены в Кельвинах, на Земле не бывает. Дело в том что данные в файлах реанализа хранятся в виде целых чисел, которые занимают гораздо меньше места чем числа с плавающей запятой. Делается это для того, чтобы файлы были меньше. Чтобы получить температуру, нужно преобразовать эти данные нехитрым способом.
Если вы посмотрите на атрибуты для переменной air, полученные при помощи ncdisp, то увидите там:
scale_factor = 0.01
Для того, чтобы получить данные в Кельвинах, нужно умножить все на scale_factor и прибавить add_offset.
Вы можете просто взять эти цифры из описания файла, либо вытянуть их при помощи следующих команд:
scale_factor = netcdf.getAtt(ncid,3,'scale_factor')
К сожалению, просто так умножить и прибавить не получится. Если вы сейчас попробуете проделать эту операцию, то получите следующую ошибку:
Дело в том, что тип переменной data сейчас int16 и на числа с плавающей запятой они умножаться не хотят. Необходимо перевести данные в тип single
На самом деле можно это сделать уже при копировании данных из переменной netCDF в переменную Матлаб. Далее переводим в Кельвины:
?, если нужно, в Цельсии
Еще немного полезной информации по этому вопросу есть в этой ветке нашего форума.
наверное должно быть
data_scaled = (data_single *scale_factor) add_off;
вместо
data_scaled = (data*scale_factor) add_off;
2RUS
Да, конечно, спасибо, Руслан, поправил!
Здравствуйте, Николай!
Невозможно скачать «подопытный» файл ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis/surface/air.sig995.2010.nc — требует логин и пароль.
В связи с тем, что правительственные учереждения США не работают, сайт с которого распространяется этот файл закрыт. Попробуйте отсюда взять http://database.rish.kyoto-u.ac.jp/arch/ncep/data/ncep.reanalysis/surface/
Спасибо! Оттуда скачалось!
Если можно, еще вопрос.
Сам я к метеорологии не имею отношения (и Матлабом тоже никогда не занимался, займусь в понедельник), но меня тесть просил помочь ему «выпотрошить» вот такой файл: http://www.cru.uea.ac.uk/cru/data/temperature/CRUTEM.4.2.0.0.anomalies.nc
С этим файлом работать в точности так же, как и с «подопытным»?
Оттуда надо вытащить температуру по координатам за последние 10 лет.
Да не за что. Выглядит вполне стандартно, так что должно быть всё то же самое. Только понятное дело не нужно к значениям ничего прибавлять и умножать, там значения сразу в кельвинах даны.
А почему решили именно матлабом открывать? Может под язык которым вы пользуетесь есть netCDF интерфейс?
Добрый день ! Хотел бы спросить у вас — а как можно найти температуру и может быть ветер по району Азовского моря ?
Тестовый файл netCDF я скачал, но было бы очень интересно получить данные по району именно Азовского моря !
Как получают эти данные по температуре слоев — это все же интерполяция, это расчетные данные ? ?ли они как то температуру измеряют именно по слоям ?
Буду благодарен за любые ссылки и комментарии по этому вопросу.
? еще вопрос — Вы перешли на Python c Matlab по какой причине ?
По причине бесплатности Python ?
Мне просто кажется что Matlab очень мощный инструмент, с него довольно сложно перейти куда то еще …
С Уважением
Сергей
sergsh
Файл глобальный, поэтому, конечно, и Азовское море там тоже есть. Однако разрешение данных будет очень грубое. На сайте NCEP есть и ветер и другие характерситики атмосферы, но если вам нужно разрешение побольше, то лучше поискать какие -либо региональные продукты.
Эти данные получают путём ассимиляции данных натурных наблюдений в модель циркуляции атмосферы. По сути можно сказать, что это динамическая инетрполяция.
Мои мысли по поводу перехода на питон я приводил в комментариях к этой записи http://koldunov.net/?p=778#comments На английском про это (со ссылками на очень подробные статьи) написано тут http://nbviewer.ipython.org/urls/raw.github.com/koldunovn/python_for_geosciences/master/00%20-%20Why%20Python.ipynb
Если вкратце, то Матлаб инструмент мощный, но не гибкий, за меняющимся миром он не поспевает.
Добрый день!
Спасибо большое за такую хорошую статью. Так как я только начал работу с данным типом файлов, очень помогла. Но возникла следующая проблема: пытаюсь считать информацию по скорости движения льда в Арктике вот с этого сайта:
ftp://ftp.ifremer.fr/ifremer/cersat/products/gridded/psi-drift/data/arctic/amsre-merged/2-daily/netcdf/2003/
?спользуя команду data = netcdf.getVar(ncid,varid); автоматически используется тип данных int16. Но вот в одном из массивов, как раз с данными о расстояниях, которые проходит лед, почти все элементы равны 32767, что означает, насколько я понимаю, выход за пределы типа данных. Хотя все элементы в этом массиве не должны превышать несколько сотен по-идее. Пробовал также читать данные с помощью команды data = netcdf.getVar(ncid,3,’single’), т.е. сразу задавать тип. Но все равно не помогает. Подскажите пожалуйста, в чем может заключаться проблема и ее возможное решение.
По вашей ссылке
ftp://ftp.ifremer.fr/ifremer/cersat/products/gridded/psi-drift/data/arctic/amsre-merged/2-daily/netcdf/2003/
открыл первый попавшийся файл 20030101-20030103.nc
и прочитал так :
ncid = netcdf.open(‘20030101-20030103.nc’,’NC_NOWRITE’);
[varname, xtype, varDimIDs, varAtts] = netcdf.inqVar(ncid,0)
varid = netcdf.inqVarID(ncid,varname)
data = netcdf.getVar(ncid,varid)
whos data
получил
data =
902880
Name Size Bytes Class Attributes
data 1×1 4 int32
то есть netcdf.getVar дает число — int32 , вроде его числом 32767 никак не переполнишь …
Стало интересно — а почему в файле 75 Кб всего одно число ?
Привет,
32767 — это отсутствующее значение, их в этих файлах должно быть дофига. Вам нужно поменять их на NaN.
Не забывайте про офсет и скейл (насколько я вижу, чтобы получить правильные значения, вам нужно просто умножить на 100)
Спасибо большое! Все получилось. Я просто думал что это ошибка какая-то, уж слишком много отсутствующих значений.
Zdrabstbyite!
Pomogite mne pojaluista, kak poluchiti dannie tsentrob deistbii atmosperii iz NCEP reanaliza c togo bremeni do nastoyashego bremeni,