Записки океанолога — обработка и визуализация данных

Задача: Сделать карту Арктики с рельефом дна, нанести на неё точки океанологических станций и нарисовать стрелочки между ними, показывающие как двигался пароход.

�?нструмент: GMT (Generic Mapping Tools).

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

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

Сайт проекта GMT. При работе всегда должна быть открыта документация, без неё никуда.

Здесь я по просьбе одного начинающего юзера расскажу о том как решить одну конкретную задачу, попутно постараюсь объяснить основы GMT. Рассказ будет довольно подробный, поскольку что такое шел скрипты (без которых GMT превращается в неповоротливое чудище) наш юзер тоже не знает.

Надеюсь что с установкой вы уже справились (там всё сравнительно безгеморойно) и мы можем переходить к вводным процедурам. Если вдруг не справились то вот тут инструкция от Алексея Колдунова для Убунту.

Первое что нам надо сделать это скачать топографию. Топография это грубо набор данных щирота\долгота\превышение на определённой сетке. Сетка может быть покрупнее или помельче. Мы выберем нечто среднее, под названием ETOPO2. Двойка означает что шаг сетки на земной поверхности — 2 минуты. Это не лучшее из доступных разрешений, но если брать более качественные то обсчитываться она будет чрезвычайно долго, а времени у нас нет ), ну и плюс для неё существует файл в нативном GMT формате .grd.

Скачать файл можно отсюда, предварительно зарегистрировавшись. Там всё оказалось нелинейно, поэтому последовательность действий следующая.

1. Регимся на сервере.
2. Переходим по ссылке
3. Внизу жмём по синей картинке с надписью Data access
4. Нажимаем Internet Download
5. Если не вошли то входим
6. Нажимаем на ETOPO2 Official Release (10800×5400 grids), September 2001
7. Скачиваем etopo2.grd.gz при помощи правой кнопки и Save as (иначе в файрфоксе, например, файл начнёт открываться в броузере, и ничего хорошего из этого не выйдет )).

Скачанный файл мы распаковываем

Получаем файл etopo2.grd. В нём содержится топография всей Земли в среднем разрешении. Нам вся Земля без надобности, мы хотим только Арктику вырезать. Для этого используем первую утилиту GMT под названием grdcut. Она вырезает определённый район из .grd файла.

• -R — определяет область данных которую мы должны вырезать. формат следующий west, east, south, and north. В нашем случае берём весь шарик по долготе, и от 70 до 90 (северного полюса) по широте. Это не вся Арктика, но чтобы сэкономить время мы пока поиграемся с этим (довольно большим) райончиком.
• -fg так толком и не понял для чего эта опция, но без неё утилита работать отказывалась.
• ./etopo2.grd — путь к нашему исходному файлу
• —G — имя выходного файла. В нашем случае arctic.grd

�?сходные данные есть, начинаем писать скрипт по отрисовке карты. Знакомые с линуксом могут пропустить следующий параграф и перейти непосредственно к скрипту. Для тех же кто не в курсе — поясняю

Создание карты в GMT это последовательность операций каждую из которых можно, конечно, набирать в командной строке, но это, учитывая их количество, меганеудобно. В таких случаях в линуксе используют шел скрипты. Грубо — это файлы, содержащие в себе последовательность команд, как если бы мы их вводили руками. Объяснять в подробностях — это долго и не ко мне 🙂 . Надеюсь что даже без познаний в шелл скриптинге всё у вас получится, я постараюсь быть предельно подробен.

В той же папке что и наш грид, создаём текстовый файл и называем его arctica.sh. Выглядеть он будет поначалу так:

А теперь давайте долго и нудно разбирать что же у нас тут понаписано 🙁

• #!/bin/sh — говорит системе какую версию шела она будет использовать. С этим проблем возникнуть не должно, потому что sh должен быть везде.
• makecpt -Cgebco > GMT_global.cpt . Утилита makecpt создаёт цветовую гамму нашей карты. В данном случае мы используем встроенную в GMT палитру gebco, которая соответствует стандартной раскраске океанов GEBCO. Для простоты сушу пока раскрашивать не будем. -С название встроенной палитры, > значит что мы перенаправляем вывод в файл, GMT_global.cpt имя файла созданной нами палитры.

grdimage arctic.grd -Js0/90/3i/70 -R-180/180/70/90 -CGMT_global.cpt > arctic.ps Тут всё самое главное и происходит.

• Утилита grdimage создаёт изображение из нашего гридированного файла.
• arctic.grd — имя гридированного файла,
• -Js0/90/3i/70 — проекция в которой будет отображаться твоя карта и её характеристики. В данном случае -Js генеральная стереографическая проекция, 0/90/ — lon0/lat0 — координаты центра проекции (центра карты), в нашем случае мы центрируемся на нулевом меридиане и северном полюсе. 3i/70 радиус (тут в дюймах) от центра до определённой широты (тут 70я).
• -CGMT_global.cpt Опция -С передаёт название файла с палитрой для карты (в нашем случае GMT_global.cpt). Знак > означает направление вывода в файл arctic.ps.

Сохранив скрипт надо дать ему права на исполнение, делается это следующим образом

теперь запустим скрипт командой

Если мы всё сделали правильно, то через какое то (довольно продолжительное) время, скрипт завершит работу и у нас образуется файлик arctic.ps

Выглядит оно довольно блекло и на карту похоже только в самых общих чертах. Давайте придадим форму этому содержанию, то есть нарисуем береговую линию и рамку карты. Дополняем наш скрипт и получаем

• В опциях grdimage добавилось -K. Это значит что в файл будет записано ещё что то.

pscoast -Js0/90/3i/70 -R-180/180/70/90 -Ba30/10 -Df -W0.1p/0/0/0 -V -O >> arctic.ps .

• Утилита pscoast отвечает за отрисовку береговой линии.
• -Js0/90/3i/70 -R-180/180/70/90 — эти опции нам уже знакомы и означают они тоже самое что и в случае grdimage.
• -Ba30/10 — это самая сложная и загадочная опция в тех утилитах GMT которые её используют 🙂 Грубо — она определяет внешний вид рамки. В данном случае из всей конструкции a30/10 я понимаю только что 30 отвечает за шаг обозначений долгот (10 видимо отвечает за широты).
• Опция -D говорит о том какое разрешение береговой линии будет использовано. Они не все ставятся по умолчанию при установке. В данном случае f означает (f)ull, то есть самое высокое разрешение. Если будет ругаться попробуйте (h)igh, (i)ntermediate, (l)ow, или на худой конец (c)rude).
• Опция -W говорит что рисовать береговую линию мы таки будем а параметры 0.1p/0/0/0 означают ширина линии/R/G/B где RGB это цветовые координаты (все нули — чёрный).
• Опция -V говорит о том что вам будут сообщаться более подробные сведения о работе утилиты (verbose mode).
• Опция -O означает что рисоваться береговая линия будет поверх уже существующей информации (батиметрии в нашем случае).
• Знаки >> означают что файл в который будет производиться запись будет не создан заново а дополнен.
• arctic.ps — имя выходного файла.

В итоге получаем гораздо более адекватный вариант, больше похожий на карту.

Простейшую карту с батиметрией мы создали. Выглядит она может и не как в атласе, но мы пока только учимся 🙂

Теперь хотелось бы нанести на карту места наших боевых подвигов, сиречь океанологических станций. Для этого нам нужен файл с их координатами (в десятых долях градуса!!!) и ещё кой какой информацией.

Создаём файл track, который выглядит следующим образом

180 90 0.2 c
31 88.33 0.2 c
31 86.67 0.2 c
31 85 0.2 c
31 83.33 0.2 c

Формат — долгота\щирота\размер значка\вид значка. По поводу размера и вида объясню позже. Видоизменяем скрипт

В pscoast появилась новая опция -K, значит будем дописывать в файл.

psxy track -Js0/90/3i/70 -R-180/180/70/90 -Gblue -S -O >> arctic.ps — рисуем станции на карте.

• psxy — утилита которая отрисовывает линии, полигоны и символы на карте.
• track — название файла с входными данными
• -Js0/90/3i/70 -R-180/180/70/90 — этих ребят мы уже знаем
• -Gblue — цвет значков, в данном случае синий, можно задавать RGB координаты.
• -S мы таки рисуем символы. В третьей и четвёртой колонке исходного файла содержится информация о размере (у нас везде 0.2) и виде (у нас c — кружок) значка. Если в файле этого нет, то можно задать эти вещи как параметры опции -S.
• -O рисуем поверх уже нарисованного.
• >> дописываем информацию в файл
• arctic.ps выходной файл.

Результат

Теперь стрелочки… Мы хотим показать куда двигался наш пароход. Для этого создаём файл vec

180 90 31 88.33
31 88.33 31 86.67
31 86.67 31 85
31 85 31 83.33

Формат следующий — долгота начальная\широта начальная\долгота конечная\широта конечная. То есть в каждой строчке координаты начала и конца нашей стрелочки. Дополняем скрипт

В первый psxy добавляется -K (не забывайте правильно ставить -K и -O, половина ошибок и стонов по поводу того что ничего не работает или работает не так — из за неправильной их расстановки).

psxy vec -Js0/90/3i/70 -R-180/180/70/90 -Gred -Svs0.01c/0.2c/0.2c -O >> arctic.ps

• psxy уже знакомая нам утилита, но теперь ей мы будем рисовать стрелочки
• vec — файл с координатами начала и конца стрелочек
• -Js0/90/3i/70 -R-180/180/70/90 — тут всё ясно
• -Gred — цвет у стрелочек будет красный
• -Svs0.01c/0.2c/0.2c Опция -Sv говорит что мы рисуем вектор, s — вместо направления и радиуса (по умолчанию) мы используем координаты начала и конца вектора, 0.01c/0.2c/0.2c — ширина вектора/длинна указателя вектора/ширина указателя вектора
• -O не забываем про эту важную штучку!
• >> дописываем в файл
• arctic.ps — выходной файл.

М.. результат тут…

.. стрелочки мне никогда не удавались )) Я тут их ещё сделал побольше чтобы было видно что это стрелочки. В общем получилось не очень, но возможность доказана 🙂

В следующий раз, когда я соберусь писать про GMT, покажу как при помощи добавления всего ещё пары строчек можно добиться вот такого результата.

Чутка на русском тут

На английском — The GMT Grdimage Primer и Introduction to Generic Mapping Tools

На вопросы которые, надеюсь, возникнут, постараюсь ответить в коментах.