GrADS и Windows. Продвинутое рисование или как сделать 2D из 4D?

Цель: Освоить новые возможности GrADS для анализа гидрометеорологической информации

Средства: GrADS

Рассмотрим GRIB-файл "model.grb". Пусть в нём содержатся данные о нескольких гидрометеорологических параметрах на нескольких вертикальных уровнях за период в несколько дней на сетке, охватывающей весь земной шар, т.е. поработаем с каким-нибудь реанализом или модельной продукцией =). Требуется получить вертикальный разрез вдоль 90 градуса восточной долготы зональной компоненты скорости ветра u. Ну и представить это в приличном для публикации виде.

Ранее мы рассматривали файл в формате netcdf. Данный формат включает в себя метаданные, т.е. данные о данных, описание данных. Другие форматы (например GRIB формат), с которыми может работать GrADS, могут не обладать таким свойством. Для того, чтобы GrADS с мог с ними работать необходимо использовать файл-описание в формате "*.ctl" (т.н. "ситээлька"). Данные были взяты с сайта GrADS Tutorial , их можно свободно скачать.

Предположим у нашего файла  "model.grb" уже есть ctl-файл (для GRIB нам понадобится ещё и "*.gmp" файл, это gribmap index file ). Выглядит он примерно так (комментарии после знака "!" авторские. В ctl-файле их нет):

dset ^model.grb   ! dset определяет какой файл будет открыт при вызове данного ctl-файла

title "Sample Model Data for lats4d Tutorial" ! Название-описание данных

undef 1e+20 ! undef определяет отсутствующие данные, пропуски

dtype grib ! тип-формат файла

index ^model.gmp ! ссылка на файл формата "*.gmp" (лежит в той же папке, где и ctl-файл). Нужен для формата GRIB

xdef 72 linear 0.000000 5.000000 ! Задаётся сетка по оси x (долгота), 72 точки линейно распределены от 0.0 c шагом в 5.0 градусов 

ydef 46 linear -90.000000 4.000000 ! Задаётся сетка по оси y (широта), 46 точки линейно распределены от -90.0 с шагом в 4.0 градуса 

zdef 7 levels ! Описание количества вертикальных уровней

1000 850 700 500 300 200 100 ! Задание значений вертикальны уровней

tdef 5 linear 0Z1jan1987 1dy ! Описание количества временных срезов, начальная дата и шаг (здесь дни)

vars 8 ! Описание количества переменных. Идёт всегда после описания широт-долгот-высот-времени. Описание выглядит так:

"имя переменной*" - количество уровней (0 если один уровень) - "другие цифры [опционально]"  - "словесное описание [опционально]"

ps        0    1,  1,  0,  0 Surface pressure [hPa]

ua        7   33,100 Eastward wind [m/s]

va        7   34,100 Northward wind [m/s]

zg        7    7,100 Geopotential height [m]

ta        7   11,100 Air Temperature [K]

hus       7   51,100 Specific humidity [kg/kg]

ts        0   11,105,  2 Surface (2m) air temperature [K]

pr        0   59,  1,  0,  0 Total precipitation rate [kg/(m^2*s)]

endvars ! Конец описания переменных

* В этом примере имена переменных отличаются от приведенных на сайте GrADS. В этом нет ничего страшного, главное - их последовательность.

Если ctl-файла нет, его можно создать самому по приведенному выше или иному примеру в любом текстовом редакторе типа Блокнота.

Теперь попробуем открыть файл  "model.grb". Для формата GRIB можно использовать команду "open" в командной строке GrADS (ga->):

open model.grb

Произойдёт ошибка! И правильно, ведь на самом деле GrADS открывает сначала не сами файлы, а их описания. Поэтому надо открывать именно ctl-файлы!

open model.ctl

Файл откроется. Теперь можно уже приступить к решению поставленной задачи. Для этого надо нарисовать вертикальный разрез (от 1000 гПа до 100 гПа) параметр ua вдоль долготы 90 в.д. за один срез по времени (пусть 3 января 1987 года).

Данные же у нас в сетке. Необходимо "вырезать" нужные нам данные. Это делается с помощью команд "set lon" и "set lat":

set lon 90 90  ! Устанавливаем  долготу с 90 в.д. по 90 в.д., т.е. вдоль одной широты

set lat  -90 90 ! Устанавливаем  широты -90 ю.ш. по 90 с.ш.

set lev 1000 100 ! Устанавливаем  вертикальные уровни от 1000 гПа до 100 гПа

set t 3 3 ! Устанавливаем  время с 3 января 1987 по 3 января 1987

Ну вот. Теперь можно рисовать ("d ua").

Добавим красоты в виде изолиний и изолиний с послойной окраской, добавим шкалу:

set gxout shaded

d ua

set gxout contour

d ua

cbarn 

Не нравится автошкала? Исправим это, сделаем её сами! Приведу пример части скрипта (N.B. всё в одиночных кавычках!), который отвечает за рисование пользовательской шкалы.

************

* Color bar*

************

'set rgb 18 0 180 240' ! команда "set rgb N r g b" позволяет определить цвет кусочка шкалы N. N - условный номер для различия получившихся цветов. r g и b - значения RGB, красного, зелёного и голубого соответственно, которые смешиваются, и в результате получается цвет N  

'set rgb 19 0 150 240' ! Т.о. задаём несколько цветов с номерами от 18 до 30. 

'set rgb 20 0 0 240'     ! Кстати, цвета от 0 до около 15 зарезервированы  под стандартные цвета. 

'set rgb 21 0 0 230'     ! Цвет 0 - бесцветный

'set rgb 22 0 0 210'

'set rgb 23 0 0 170'

'set rgb 24 0 0 130'

'set rgb 25 0 0 90'

'set rgb 26 90 0 0'

'set rgb 27 130 0 0'

'set rgb 28 170 0 0'

'set rgb 29 210 0 0'

'set rgb 30 250 0 0'

'set clevs -10 -5 0 5 11 15 20 25 30 35 40' ! Команда "set clevs" задаёт цену деления (могут быть неодинаковыми, но будут рисоваться слева направо) будущей шкалы. В GrADS шкалы, которые получаются в результате вызова cbarn и cbar, имеют на каждом конце значения "<" и ">" самого маленького и самого большого из установленных значений. Поэтому на P значений clevs приходится P+2 значений ccols. 

'set ccols 18 19 20 21 22 23 24 25 26 27 28 29 30'

Команда "set ccols" задаёт последовательность цветов по условным номерам (подготовленным заранее или стандартным). Промежуткам между значениями clevs (начинается с "меньше самого малого") соответствуют значения ccols.

Снова рисуем ("d ua") поле. Теперь картинка хороша, но хочется ещё подписать оси и название рисунка. Тут всё просто, используем команды "draw title", "draw xlab" и "draw ylab":

'draw title Zonal average U-wind, m/s' ! Рисуем заголовок

'draw xlab Latitude, grad' ! Рисуем подпись оси x

'draw ylab Vertical Levels , hPa' ! Рисуем подпись оси y

А теперь усложним задачу. Пусть мы хотим нарисовать вертикальный разрез ОСРЕДНЕННЫХ по времени и долготам величины ua! Тут нам понадобится команда осреднения "ave". Мы сначала осредним значения ua по долготам, а потом по времени:

ga-> define U2 = ave(ave(ua,lon=90,lon=360),t=1,t=5)

Это эквивалентно последовательности команд

ga-> define U2a = ave(ua,lon=90,lon=360)

ga-> define U2   = ave(U2a,t=1,t=5)

В GrADS очень неплохой скриптовый язык, в нём есть даже операторы ветвления и цикла. Так вот, с помощью команды "define" можно определить новую переменную и записать в неё массив данных. Теперь мы просто рисуем переменную U2:

d U2

Вся получившаяся красота занимает много места, набирать столько команд каждый раз неудобно. Запишем всё это в скрипт "Vert_profile.gs" ($N - примечания, в листинге скрипта их нет):

********************** $1

* Start grads script *

**********************

'reinit'

'open model.ctl'

*************

* Set block *

*************

lev1=1000 $2

lev2=100

lat1=-90

lat2=90

lon1=90

lon2=360

U=ua $3

'set lon 'lon1' 'lon1'' $4

'set lat 'lat1' 'lat2''

'set lev 'lev1' 'lev2''

------------------------------------------------------------

************

* Color bar*

************

'set rgb 18 0 180 240'

'set rgb 19 0 150 240'

'set rgb 20 0 0 240'

'set rgb 21 0 0 230'

'set rgb 22 0 0 210'

'set rgb 23 0 0 170'

'set rgb 24 0 0 130'

'set rgb 25 0 0 90'

'set rgb 26 90 0 0'

'set rgb 27 130 0 0'

'set rgb 28 170 0 0'

'set rgb 29 210 0 0'

'set rgb 30 250 0 0'

'set clevs -10 -5 0 5 11 15 20 25 30 35 40'

'set ccols 18 19 20 21 22 23 24 25 26 27 28 29 30'

------------------------------------------------------------

****************************

* Averaging & Ploting block*

****************************

'define Uave = ave(ave('U',lon='lon1',lon='lon2'),t=1,t=5)'

****************

* Ploting block*

****************

'set grads off'

'set csmooth on'

'set mpdset hires'

'set gxout shaded' 

'd Uave'

'set gxout contour'

'set cthick 6'  $5

'd Uave'

'cbarn 1 0 ' * 1 - нормальный размер, 0 - горизонтально

'draw title Zonal average U-wind, m/s'

'draw xlab Latitude, grad'

'draw ylab Vertical Levels , hPa'

'printim vert_profile.gif gif x800 y600 white'

Получившийся результат:

Здесь я использовал ещё некоторые команды (они помечены цветом-знаком $):

1) строки, которые не могут быть выполнены и не приводят к явным ошибкам, игнорируются. Их можно использовать как комментарии (я использую знаки * и # )

2) вводится переменная lev1, которой присваивается значение 100. Обратите внимание, здесь нет кавычек!

3) вводится переменная U, которой передается переменная-массив ua

4) вместо числовых констант при определении границ области (set lon-lat) командам передаются введенные выше переменные. Обратите внимание на кавычки! Переменные должны быть в кавычках, когда вы используете их в командах (сами команды в целом также в кавычках).

5) команда "set cthick" устанавливает толщину изолиний (6 сильно жирнее, чем 1 и 4)

В примерах GrADS показаны также и другие возможности этой системы с использованием тех же данных.

P.S. Часть скрипта (помечена ---), которая "рисует" пользовательскую шкалу, можно сохранить как отдельный скрипт "my_scale.gs" и вызывать его командой run. Получится матрёшка. =)