Функции NetCDF в Мatlab 7.9 (R2009b)

Большую часть времени у океанолога занимает общение с различными программами и написание небольших программ для обработки и визуализации данных. Обсуждению этого и посвящён данный раздел.

Модераторы: Magik, koldunov aleksey, kostek

Функции NetCDF в Мatlab 7.9 (R2009b)

Сообщение Dima » Ср окт 07, 2009 12:41 pm

Помогите, пожалуйста, разобраться, как можно прочесть файл данных реанализа температуры (ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis.dailyavgs/surface/air.sig995.2008.nc) при помощи встроенных функций обработки netcdf, которые поддерживает новая версия Matlaba? В частности необходимо построить распред. температуры для любого заданного дня года.

Код: Выделить всё
% Open file.
ncid = netcdf.open('air.sig995.2008.nc','NC_NOWRITE');

% Get the name of the first variable.
n=2
[varname, xtype, varDimIDs, varAtts] = netcdf.inqVar(ncid,n);

% Get variable ID of the first variable, given its name.
varid = netcdf.inqVarID(ncid,varname);

% Get the value of the first variable, given its ID.
data = netcdf.getVar(ncid,varid)


Проблемы в следующем:после выполнения выше написанных функций, массивы дат и температур выводится в каком-то непонятном числовом виде непригодном для дальнейшей обработки, а широта и долгота -правильно.
А нужно чтобы выводилась матрица распред. температуры для любого заданного дня года.
Аватара пользователя
Dima
Завсегдатай
 
Сообщения: 103
Зарегистрирован: Пт окт 02, 2009 9:05 am
Откуда: Москва

Сообщение Magik » Чт окт 08, 2009 12:54 pm

Проблема ваша - всегдашняя проблема людей которые начинают работать с реанализом ) Если вы посмотрите на дамп заголовка файла, то увидите следующее (в винде вроде тоже должно работать из командной строки)

Код: Выделить всё
ncdump -h air.sig995.2008.nc
netcdf air.sig995.2008 {
dimensions:
        lon = 144 ;
        lat = 73 ;
        time = UNLIMITED ; // (366 currently)
variables:
        float lat(lat) ;
                lat:units = "degrees_north" ;
                lat:actual_range = 90.f, -90.f ;
                lat:long_name = "Latitude" ;
        float lon(lon) ;
                lon:units = "degrees_east" ;
                lon:long_name = "Longitude" ;
                lon:actual_range = 0.f, 357.5f ;
        double time(time) ;
                time:units = "hours since 1-1-1 00:00:0.0" ;
                time:long_name = "Time" ;
                time:actual_range = 17593032., 17601792. ;
                time:delta_t = "0000-00-01 00:00:00" ;
                time:avg_period = "0000-00-01 00:00:00" ;
        short air(time, lat, lon) ;
                air:long_name = "mean Daily Air temperature at sigma level 995" ;
                air:valid_range = 185.16f, 331.16f ;
                air:actual_range = 189.84f, 316.47f ;
                air:units = "degK" ;
                air:add_offset = 512.81f ;
                air:scale_factor = 0.01f ;
                air:missing_value = 32766s ;
                air:precision = 2s ;
                air:least_significant_digit = 1s ;
                air:GRIB_id = 11s ;
                air:GRIB_name = "TMP" ;
                air:var_desc = "Air temperature\n",
                        "A" ;
                air:dataset = "NCEP Reanalysis Daily Averages" ;
                air:level_desc = "Surface\n",
                        "0" ;
                air:statistic = "Mean\n",
                        "M" ;
                air:parent_stat = "Individual Obs\n",
                        "I" ;

// global attributes:
                :Conventions = "COARDS" ;
                :title = "mean daily NMC reanalysis (2008)" ;
                :history = "created 2007/12 by Hoop (netCDF2.3)" ;
                :description = "Data is from NMC initialized reanalysis\n",
                        "(4x/day).  These are the 0.9950 sigma level values." ;
                :platform = "Model" ;
}


важными для вашего случая являются строчки
Код: Выделить всё
air:add_offset = 512.81f ;
air:scale_factor = 0.01f ;


Дело в том что для уменьшения объема файлов данные по температуре храняться в формате short и для того чтобы перевести их в удобоваримый вид нужно умножить на scale_factor (в вашем случае 0.01) и добавить офсет (в вашем случае 512.81). Для полного счастья можно еще перевести кельвины в цельсии произведя вычитание 273.15 :) .

Подробнее о работе с реанализом можно почитать например тут
http://koldunov.net/?p=310
возможно вы захотите использовать что то более приспособленное для работы с NetCDF чем матлаб )
Аватара пользователя
Magik
Site Admin
 
Сообщения: 2133
Зарегистрирован: Чт авг 24, 2006 9:21 pm
Откуда: Питер/Гамбург

Сообщение koldunov aleksey » Чт окт 08, 2009 2:16 pm

У меня не такая новая версия Матлаба и я не в курсе их новых функий обработки netcdf. Я открыл Ваш файл командой:
Код: Выделить всё
ncload ('air.sig995.2008.nc')

для удобства использования и построения изображения переопределил матрицу к виду (широта, долгота, время) командой:
Код: Выделить всё
air = permute (air, [2 3 1]);

для простейшего построения (без привязки к картографической основе) можно использовать команду:
Код: Выделить всё
imagesc (air(:,:,5))

(где "5" - срок наблюдения, т.е. его порядковый номер в матрице).
Перевести температуру к человеческим цельсиям можно как описано выше
Аватара пользователя
koldunov aleksey
Коренной Житель
 
Сообщения: 364
Зарегистрирован: Чт авг 24, 2006 10:22 pm
Откуда: Санкт-Петербург

Сообщение Dima » Чт окт 08, 2009 2:57 pm

Magik писал(а):Дело в том что для уменьшения объема файлов данные по температуре храняться в формате short и для того чтобы перевести их в удобоваримый вид нужно умножить на scale_factor (в вашем случае 0.01) и добавить офсет (в вашем случае 512.81). Для полного счастья можно еще перевести кельвины в цельсии произведя вычитание 273.15 :) .


Спасибо! Хорошо, с температурой, как переводить понятно. Но, неужели сам Матлаб не поддерживает формат short и это все надо делать вручную? И что делать со временем, его же тоже надо привести к нормальному виду?
Аватара пользователя
Dima
Завсегдатай
 
Сообщения: 103
Зарегистрирован: Пт окт 02, 2009 9:05 am
Откуда: Москва

Сообщение Magik » Чт окт 08, 2009 3:34 pm

Dima писал(а):Спасибо! Хорошо, с температурой, как переводить понятно. Но, неужели сам Матлаб не поддерживает формат short и это все надо делать вручную? И что делать со временем, его же тоже надо привести к нормальному виду?


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

Magik писал(а):Обратите внимание на атрибуты add_offset и scale_factor. Они означают что для того чтобы получить действительно значения в Кельвинах, нам нужно с данными проделать следующую операцию:

Код: Выделить всё
surf_air_float = surf_air*scale_factor + add_offset


Это обычный приём применяемый для того чтобы уменьшить размер файла. Возможно вы заметили что тип переменной air в исходном файле, а теперь и нашей переменной surf_air был short, что значит её размер всего 16 бит и она может содержать значения +/- (32767). Для того чтобы хранить переменные типа float, для каждой из них требуется 32 бита. Может для небольших файлов это и не так значимо, но когда счёт идёт на террабайты, приходится экономить каждый бит. Файл довольно сильно уменьшается, при этом переход от short к float происходит чрезвычайно просто.


Таким образом это просто трюк который позволяет уместить больше информации в меньший объем файла. Только некоторые программы для отображения данных в формате NetCDF производят эту трансформацию автоматически.

Что касается времени, то да его тоже нужно приводить к нормальному виду, сейчас у вас значения в единицах:

Код: Выделить всё
time:units = "hours since 1-1-1 00:00:0.0"


то есть в часах с 1-1-1, нужно это дело переводить в нормальное время. Поищите на matalbcentral.com наверняка уже кто то сделал подобный конвертер.
Последний раз редактировалось Magik Чт окт 08, 2009 3:39 pm, всего редактировалось 1 раз.
Аватара пользователя
Magik
Site Admin
 
Сообщения: 2133
Зарегистрирован: Чт авг 24, 2006 9:21 pm
Откуда: Питер/Гамбург

Сообщение koldunov aleksey » Чт окт 08, 2009 3:38 pm

Dima писал(а): И что делать со временем, его же тоже надо привести к нормальному виду?

Про форматы времени в матлаб написано здесь:
http://www.oceanographers.ru/forum/viewtopic.php?t=777
прочитайте, если появятся вопросы - задавайте.
Аватара пользователя
koldunov aleksey
Коренной Житель
 
Сообщения: 364
Зарегистрирован: Чт авг 24, 2006 10:22 pm
Откуда: Санкт-Петербург

Сообщение Dima » Чт окт 08, 2009 5:48 pm

koldunov aleksey писал(а):Про форматы времени в матлаб написано здесь:
http://www.oceanographers.ru/forum/viewtopic.php?t=777
прочитайте, если появятся вопросы - задавайте.

Дело в том, что функция Матлаба
Код: Выделить всё
datestr
не верно конвертирует эти числовые значения для любого взятого из таблицы формата дат, т.ч. видимо даты тоже надо будет как-то вручную приводить к нормальному виду из числового вида( 17593032). Как расшифровать этот числовой формат дат netcdf?
Аватара пользователя
Dima
Завсегдатай
 
Сообщения: 103
Зарегистрирован: Пт окт 02, 2009 9:05 am
Откуда: Москва

Сообщение koldunov aleksey » Чт окт 08, 2009 7:26 pm

Dima писал(а):Дело в том, что функция Матлаба
Код: Выделить всё
datestr
не верно конвертирует эти числовые значения для любого взятого из таблицы формата дат, т.ч. видимо даты тоже надо будет как-то вручную приводить к нормальному виду из числового вида( 17593032). Как расшифровать этот числовой формат дат netcdf?

Ну Вы же знаете на какую дату у Вас первое значение (Вы же скачивали для какого-то периода). Узнайте для этого момента времени дату в матлабовском формате, вычислите разницу, и эту разницу отнимите от каждого вашего значения даты. Получится время в Матлабовском формате.
Аватара пользователя
koldunov aleksey
Коренной Житель
 
Сообщения: 364
Зарегистрирован: Чт авг 24, 2006 10:22 pm
Откуда: Санкт-Петербург

Сообщение Dima » Пт окт 09, 2009 7:29 am

koldunov aleksey писал(а):Ну Вы же знаете на какую дату у Вас первое значение (Вы же скачивали для какого-то периода). Узнайте для этого момента времени дату в матлабовском формате, вычислите разницу, и эту разницу отнимите от каждого вашего значения даты. Получится время в Матлабовском формате.


Спасибо! С датой все ясно, чтобы привести ее к нормальному виду, нужно сделать следующее преобразование:
Код: Выделить всё
real_time=datestr(time_netcdf/24+366)


Теперь, чтобы получить матрицу температур, например , для 1 января 2008 года делаем так:
Код: Выделить всё
temp = netcdf.getVar(ncid,varid,[0 0 1], [144 73 1], 'short')



Еще вопрос, может ли Матлаб открывать и брать данные из файлов nc, лежащих в сети?
Как наиболее простым образом преобразовать матрицу размерности 144 на 73 в 180 на 360?
Аватара пользователя
Dima
Завсегдатай
 
Сообщения: 103
Зарегистрирован: Пт окт 02, 2009 9:05 am
Откуда: Москва

Сообщение Magik » Пт окт 09, 2009 2:39 pm

Dima писал(а):Еще вопрос, может ли Матлаб открывать и брать данные из файлов nc, лежащих в сети?
Как наиболее простым образом преобразовать матрицу размерности 144 на 73 в 180 на 360?


По первому вопросу можно почитать тут
http://www.oceanographers.ru/forum/viewtopic.php?t=216

По второму - если вас интересует интерполяция на одноградусную сетку, то попробуйте посмотреть в сторону функции interp2
Аватара пользователя
Magik
Site Admin
 
Сообщения: 2133
Зарегистрирован: Чт авг 24, 2006 9:21 pm
Откуда: Питер/Гамбург


Вернуться в Программы и программирование