Записки океанолога — обработка и визуализация данных
Данные: значения на сетке в формате netCDF. Сетка покрывает всю Землю.
Задача: проинтерполировать эти данные на другую сетку (T63), чтобы их можно было сравнить с другой моделью.
�?нструменты: cdo и GMT для отображения.
cdo (climate data operators) это коллекция строковых операторов для обработки и анализа данных климатических моделей (и вообще любых данных в netCDF формате), распространяемая по лицензии GPL. Поддерживаемые форматы данных GRIB, netCDF, SERVICE, EXTRA и IEG. Доступно более 350 операторов.
cdo могут делать чудеса, конечно, в том что касается работы с netCDF, но то, что они могут так просто интерполировать с сетки на сетку добило меня окончательно ) Возможно, что ещё кого-то это также восхитит 🙂
Официальный сайт cdo вот здесь вот . 🙂
Ликбез
Про формат netCDF на русском можно прочитать тут
Начальные сведения о cdo на русском (инсталяция , простейшие операторы) вот тут �?сходные данные у нас на очень хитропопой сетке. Выглядит она следующим образом (смотрим на землю сверху, по центру — северный полюс):
Как видим дырка у неё не над северным полюсом, а смещена в район Гренландии. Это сделано для того, чтобы избежать схождения меридиан над предполагаемой российской территорией и забыть о разного рода трудностях с этим связанных. Поскольку ни над ни под Гренландией океана нет, решили ей пожертвовать в научных целях, несмотря на неполиткорректность этого шага.
Стандартный же вариант это, так называемая, Гауссовская сетка. Выглядит она так:
Как видим тут свет сходится клином на северном полюсе и даже суперкомпьютерам приходится тяжко при обработке этой кучи маленьких ячеек, в которых должны интегрироваться те же уравнения, что и в больших. Эту сетку в основном используют в атмосферных моделях, потому что воздух, к сожалению, есть везде, даже над Гренландией. Но есть и модели океана на такой сетке.
Так вот, чтобы сравнить данные одной модели с данными другой по пространству, мы должны с одной сетки интерполировать значения на другую.
Данные это концентрация льда и в оригинале (на хитропопой сетке) они выглядят следующим образом:
Видно, что Гренландия стыдливо прикрывает белое пятно без данных. Чтобы не заморачиваться, легенду тут не привожу, в принципе не так уж важно, что это за данные — просто распределение характеристики в пространстве.
Перевод на гауссовскую сетку при помощи cdo будет выглядеть следующим образом
- cdo — собственно вызывает cdo
- remapbic — интерполирует данные на другую сетку при помощи бикубической интерполяции (есть ещё три других варианта, описанных ниже)
- t63grid — название сетки на которую интерполируем. T63 стандартная и «зашита» в cdo (есть ещё масса вариантов описанных ниже)
- input.nc — входной файл
- output.nc — выходной файл
Всё. Гениально и просто ) В итоге получаем
Не очень красиво, согласен… Дырки чёрные какие то появились (районы без данных)… Можем попробовать другой метод интерполяции.
�?х в cdo доступно четыре штуки:
- remapbil — билинейная интерполяция
- remapbic — уже знакомая нам бикубическая интерполяция
- remapcon — conservative remapping (затрудняюсь перевести, сохраняет общее количество параметра, насколько я понимаю)
- remapdis — средневзвешенное осреднение ближайших четырёх точек.
Результат для одних и тех же данных будет получаться немного разный:
remapbil
remapdis
remapcon
Конечно, интерполировать можно не только на гауссовскую сетку. В качестве целевой сетки (мы интерполировали на гауссовсукую T63 и задали её как t63grid) могут выступать:
- rNXxNY — Глобальная регулярная сетка. В отличие от Гауссовской, у которой широты распределены с различным шагом, в регулярной распределение их равномерное, задаётся количество ячеек сетки по долготе NX и широте NY
- tRESgrid Глобальная гауссовская сетка. RES — триангуляционное разрешение
- gmeNI — глобальная гексогональная сетка. NI — количество интервалов на основной стороне треугольника (что бы это не значило )))
- file — Если грид у вас не стандартный (как тот, что у нас, с дыркой под Гренландией) вы можете описать его при помощи файла в определённом формате. Форматы описаны в документации к cdo.
Пример гексагональной сетки 🙂
Вот так вот, без выгибонов мы интерполируем данные туда сюда )
На последок скрипты при помощи которых я отрисовывал сетку и данные (рисовалка GMT). Там используется недокументированная команда cdo — outputbounds, которая переводит данные netCDF в ASCII файл, формата понимаемого GMT.
set proj = "s0/90/3i/70" set regi = "-180/180/70/90" set setup = "-J $proj -R $regi "
set epsname = T63_grid.ps
cdo outputbounds S09_air2m_mm_1980- 1999 _mon01_mean.nc > ASCII_FILE.txt
psxy ASCII_FILE.txt $setup -Y3 -K -M -L -P > $epsname
pscoast $setup -Ba30 / 10 -W0.1p / 0 / 0 / 0 -O -P >> $epsname
set proj = "s0/90/3i/70" set regi = "-180/180/70/90" set setup = "-J $proj -R $regi "
makecpt -Crainbow -T0 / 1 / 0.1 > color1.cpt
set palet = color1.cpt set epsname = remapcdis.ps set ncname = remapdis.nc
cdo outputbounds $ncname > ASCII_FILE.txt
psxy ASCII_FILE.txt $setup -C $ -Y3 -K -M -L -P > $epsname
pscoast $setup -Ba30 / 10 -G150 -O -P >> $epsname
6 комментариев
а какую с какой моделью ты собрался сравнивать?
2Sca Это просто для примера ) Вообще я сравниваю спутниковые данные с моделью. Спутниковые интерполируются немного геморойней и сложнее было бы объяснить принцип. )
Виктор, а Вы так на вскидку не скажете — возможно ли в CDO интерполировать данные скажем из сетки в географических координатах в сетку в геомагнитных? Если такая возможность есть (пускай даже ручками как-то задавать сдвиг оси и т.д.), то есть смысл тоже использовать CDO как инструмент!
2Konstantin Виктор это кто? 🙂 Я не очень себе представляю как выглядит сетка геомагнитных координат, но если ее можно описать в виде — координаты центра\координаты границ, то почему бы и нет.