Визуализация данных при помощи PyNGL
Задача: взять netCDF файл содержащий две компоненты скорости течений (плюс, естественно, координаты) и изобразить вот это:
Инструменты: PyNGL и PyNIO
Смотрю я на то чем пользуется океанологическая общественность вокруг и замечаю что всё больше и больше народу тянуться к питону. У матлаба и ему подобных есть много преимуществ, но не меньше недостатков, главный из которых их небесплатность. Питон же бесплатен и к тому-же является серьёзным языком программирования на котором пишут не только научные работники.
Ну и начать хотелось с чего то осязаемого. Ничего более осязаемого чем визуализация быть не может, поэтому для бесчеловечных экспериментов был выбран свободно распространяемый модуль PyNGL (пингль), рассчитанный на людей занимающихся геонауками. Преимущества пингля над тем же в том что данные можно считывать из netCDF (и кучи других форматов) непосредственно и использовать для обработки этих данных весь арсенал SciPy и питона как такового.
Сразу оговорюсь что в питоне я мегановичок (вчера прочитал по нему первую статью) так что если какие то вещи покажутся очевидными, смело пропускайте. Поехали.
Для начала опишем инструменты:
это модуль для питона использующийся для визуализации научных данных с упором на высококачественную 2D визуализацию.
PyNIO это тоже питоновски модуль, поставляемый вместе с PyNGL, и использующийся для чтения и записи файлов в нескольких форматах имеющих широкое распространение в геонауках. При помощи модуля можно работать с netCDF, GRIB1, HDF 4, HDFEOS 2, CCM history files, GRIB2, netCDF 4.
Вывод данных производится по выбору в PostScript, PDF, X11 screen или NCGM (этот вам врятли когда понадобится).
Происходит пингль от NCAR Command Language(NCL) который разрабатывается в NCAR и которому, по слухам, всех работников NCEP и NCAR заставляют учиться (Маша, привет :)).
Установка PyNGL
Для начала вам нужен Python. Как его ставить объяснять не буду :). Затем вам необходимо озаботиться установкой модуля - предназначенный для научного программирования, который потянет NumPy в качестве зависимости.
Делается это в Убунте следующим образом:
Если же вы всё таки хотите только NumPy то
Дальше начинается небольшая канитель. PyNGL/PyNIO поставляются в виде бинарников для различных систем (линукс, виндоуз, макось и солярис) и загружаются с сайта . Для того чтобы загрузить нужно зарегистрироваться. Для этого ползём на сайт ESG и смотрим в левую часть под "ESG News". Там есть кнопочка запросить аккаунт. В поле "Statement of work" укажите что вы "intrested in using PyNGL/PyNIO". После этого придёт автоматическое подтверждение того что заявка получена. Заявка рассматривается вручную админом. У меня ушло на это 2 часа. Могут рассматривать до 3 рабочих дней :).
Скачиваете бинарник для своей системы и устанавливаете его:
tar -xf PyNGL-1.2.0.Linux-i686-gcc4-2.5.tar
cd PyNGL-1.2.0
sudo python setup.py install
Если вы работаете под Ubuntu то новую версию 1.3b ставить я пока не советую. выбираейте стабильную 1.2.
Для проверки того всё ли нормально можно выполнить
Эта команда выполняет автоматом скрипты уже написанные для нас добрыми людьми. Список примеров доступен .
Ну и последнее. Для того чтобы работать с программой я буду использовать . Это такая надстройка, что ли, над питоновским шелом которая делает жизнь немного приятнее. Вы можете поставить его чтобы в точности следовать происходящему в этой статье, но для исполнения скриптов вам он не потребуется.
Установка в убунту, как не трудно догадаться:
Я буду следовать примеру расположенному , и продвинутая публика возможно предпочтёт туда перейти сразу. Те же кто, как и я, в питоне и пингле не бум бум, давайте разбираться вместе.
Начинаем писать скрипт
По сто раз набивать одни и те же команды, конечно не с руки, так что сразу начнём писать скрипт. Создаём в нашем любимом редакторе с подсветкой (в моём случае Kate) текстовый файл и называем его vector.py. Вставляем следующий код:
# Detiled discription in Russian you can find at http://koldunov.net/?p=148
# Script is based on http://www.pyngl.ucar.edu/Examples/Scripts/vector_pop.py
import Ngl
import Nio
import numpy
#
# Open a workstation.
#
rlist = Ngl.Resources()
rlist.wkColorMap = ["White","Black","Tan1","SkyBlue","Red"]
wks_type = "ps"
wks = Ngl.open_wks(wks_type,"vector",rlist)
где:
- import - импортирует модули необходимые для работы скрипта.
- всё остальное - не должно вас занимать. Эти команды готовят "почву" для нашей работы. Единственно что тут интересного пока это название выходного файла - vector.
После этого заходим в ipython и запускаем скрипт:
(возможно перед run вам необходимо будет поставить знак %). Если вы не работает в ipython то можно просто набрать
Дальше все команды будут относиться к ipython.
Чтобы посмотреть на то что мы натворили (на наши переменные) нужно в командной строке ipython набрать:
получим:
То етсь мы загрузили три модуля и насоздавали переменных, углубляться в которые нам сейчас не обязательно.
Загрузка данных
Теперь нужно дать скрипту данные. В наличии у нас файл netCDF с двумя компонентами скорости и координатами, одна из которых время. . Надеюсь что библиотеки netCDF и сопутствующие программы у вас уже установлены и вы можете посмотреть на него:
Получиться должно вот что:
Для тех кто хочет поподробней понять что тут да как отсылаю к небольшому . Основная же информация тут следующая - у нас имеется 6 "переменных":
- lon - долгота, матрица имеющая размеры y на x, в нашем случае 220x254.
- lat - широта, всё то же самое что и для долготы
- lev - уровень (глубина), эта переменная зависит от измерения lev (но нам на это наплевать сейчас)
- time - время, зависит от времени (даже не думайте об этом ))
- SICUO - в этой переменной хранятся u составляющие скорости. Размеры у неё (time, lev, y,x). Это важно )
- SICVE - всё то же самое только для составляющей v.
Информацию к сведению приняли, теперь добавляем в скрипт строчки
# Open the netCDF file.
#
file = Nio.open_file("./S09_velos_09.nc")
#
# Get the u/v and lat/lon variables.
#
urot = file.variables["SICUO"]
vrot = file.variables["SICVE"]
lat2d = file.variables["lat"]
lon2d = file.variables["lon"]
где
- file = Nio.open_file("./S09_velos_09.nc") - открывает наш файл, он должен находиться в той же директории что и скрипт.
- urot = file.variables["SICUO"] - и все остальные ассоциируют переменные (urot, vrot и т.д.) с переменными в netCDF файле.
запускаем и выполняем whos. Получаем
Так как язык у нас объектно ориентированный то переменные являются объектами. Например переменная (объект) file имеет тип NioFile. Это значит что объект данного типа имеет определённые атрибуты, описаные . Мы уже использовали атрибут variables для того чтобы проасоциировать переменные (urot, vrot и т.д.) с переменными в netCDF файле. Мы также можем посмотреть какие в принципе есть переменные в файле (более лёгкий способ чем делать ncdump -h), а также какие есть в нём измерения и какие у этих измерений размеры:
Посмотрим теперь на переменную urot, тип которой - NioVariable. У этого объекта имеются следующие атрибуты:
- rank - количество измерений в переменной
- shape - количество элементов в каждом измерении
- dimensions - название измерений
- __dict__ - атрибуты асоциированные с переменной в файле
В применении к urot:
Самое главное для нас тут то что urot (как и vrot) имеют 4 измерения, а для того чтобы рисовать наши векторы нам нужно всего два. К тому же данные, насколько я понимаю, должны быть типа ndarray, поэтому кое что сделаем и с нашими переменными lat2d и lon2d, у которых и так два измерения, но тип данных другой. Добавляем в скрипт:
v = vrot[0,0,:,:]
lon = lon2d[:]
lat = lat2d[:]
- u = urot[0,0,:,:] - помните что измерения у urot были ('time', 'lev', 'y', 'x'). Первые два нас не интересуют, поэтому на их место мы ставим нули. Знак двоеточия как и в матлабе означает "все элементы". В итоге данная запись означает - создать переменную u со всеми элементами из двух последних измерений ('y', 'x') переменной urot.
- lon = lon2d[:] - в новую переменную импортируется весь 2D массив (оба измерения).
Исполняем, и вводим команду whos:
Самое сложное - импортирование данных, позади. Теперь осталась одна развлекуха с рисованием.
Визуализация данных
Перво наперво мы должны создать переменную с типом Ngl.Resources в атрибутах которой будут содержаться все сведения касающиеся нашего рисунка. Подробнее о ресурсах . Добавляем в скрипт:
## Now set up a vector resource variable.
##
vcres = Ngl.Resources()
##
## Set coordinate arrays for data. This is necessary in order to overlay
## the vectors on a map.
##
vcres.vfXArray = lon
vcres.vfYArray = lat
vcres.mpProjection = "Orthographic"
vcres.mpFillOn = True
vcres.mpLandFillColor = "Tan1"
vcres.mpOceanFillColor = "SkyBlue"
vcres.mpInlandWaterFillColor = "SkyBlue"
vcres.mpLimitMode = "LatLon"
vcres.mpCenterLonF = -80
vcres.mpCenterLatF = 55
vcres.mpMinLatF = 60
vcres.mpDataBaseVersion = "MediumRes"
vcres.vcRefAnnoOrthogonalPosF = -0.18 # Move ref anno up into plot.
vcres.tiMainString = "Currents"
plot = Ngl.vector_map(wks,u,v,vcres)
- vcres = Ngl.Resources() - переменную с типом Ngl.Resources
- vcres.vfXArray = lon - задаём значение ресурса в котором содержатся значения долгот
- vcres.vfYArray = lat - то же для широт
- vcres.mpProjection - проекция, понятно дело
- vcres.mpFillOn - заливать или не заливать области на карте (сушу и океан)
- vcres.mpLandFillColor - цвет заливки земли
- vcres.mpOceanFillColor - цвет заливки океана
- vcres.mpInlandWaterFillColor - цвет для речек и озёр
- vcres.mpLimitMode - определяет как будет выбираться видимая область карты. В нашем случае мы зададим её при помощи широт и долгот.
- vcres.mpCenterLonF - долгота центра проекции
- vcres.mpCenterLatF - широта центра проекции
- vcres.mpMinLatF - минимальная широта карты
- vcres.mpDataBaseVersion - качество отображения берегов и политических границ.
- vcres.vcRefAnnoOrthogonalPosF - контролирует положение вашего референс вектора (легенды грубо говоря)
- vcres.tiMainString - название вашего графика
Описание всех ресурсов находится . Первые, маленькие, буквы говорят о том какой это тип ресурса. Например vc - управляет свойствами вектора, mp - свойствами карты, ti - заголовками.
Наконец, мы рисуем нашу картинку командой plot = Ngl.vector_map(wks,u,v,vcres), где
- wks - это наша рабочая станция, параметры которой задавались в самом начале (вывод в .ps, используемые цвета и проч)
- u - u компонента течения, которую мы так нудно выковыривали
- v - v компонента
- vcres - удивительная переменная в которой храняться все параметры (значения ресурсов) относящиеся к нашей карте.
Запускаем скрипт (питон может слегка поругаться, игнорируем) и в итоге должны получить в директории со скриптом файл vector.ps:
Цвета приятные но информации не много 
Это модельные данные на специальной сетке полюс у которой расположен над (или под) Гренландией чтобы избежать проблем со схождением меридиан у географического полюса. Отчётливо видно что точек гораздо больше около модельного полюса. Но разобрать ничего не возможно, слишком много точечек.
Попробуем исправить положение отлисовывая только каждый третий вектор. Добавляем в скрипт:
## Now we're going to thin the vectors by drawing only every third one.
##
vcres.vfXArray = lon[::3,::3]
vcres.vfYArray = lat[::3,::3]
vcres.tiMainString = "Striding the grid [::3,::3]"
plot = Ngl.vector_map(wks,u[::3,::3],v[::3,::3],vcres)
Тут по идее должно быть всё понятно. Просо берём каждое третье значение по x и по y.
В вашем .ps файле появится вторая страница, а на ней вы должны увидеть что то типа:
Такого же эффекта можно добиться контролируя расстояние между векторами. Также мы изменим референс вектор. Добавляем в скрипт:
## This time thin the vectors by setting the vcMinDistanceF resource.
## We change the leight of the vectors as well
##
vcres.vfXArray = lon
vcres.vfYArray = lat
vcres.vcRefMagnitudeF = 0.3
vcres.vcMinDistanceF = 0.015
vcres.vcRefLengthF = 0.08
vcres.tiMainString = "vcMinDistanceF = 0.015"
plot = Ngl.vector_map(wks,u,v,vcres)
Возвращаем vcres.vfXArray и vcres.vfYArray изначальные значения (все точки без исключения).
Векторы значения которых равны vcRefMagnitudeF будут сделаны длинны vcRefLengthF. Короче говоря это способ контролировать длину ваших векторов и их правильное отображение. Если значение vcRefMagnitudeF будет очень большим то длинна ваших векторов будет одинаковой и не будет отражать силы течения, векторы просто покажут направления. Подбирать vcRefMagnitudeF нужно исходя из значений ваших данных (cdo info команда для тех кто использует ).
Получили:
Теперь векторы отражают силу течения, но всё равно какие-то маловатые. Давайте сделаем их побольще. Добавляем в скрипт:
## Make the vectors a little longer.
##
vcres.vcMinFracLengthF = 0.05
vcres.tiMainString = "vcMinFracLenF = 0.05"
plot = Ngl.vector_map(wks,u,v,vcres)
Значение ресурса vcMinFracLengthF задаётся в частях от вашего референс вектора vcRefMagnitudeF и означает длину самого маленького вектора. Длины всех остальных векторов распределяются между длинной этого, самого маленького и длинной самого большого - нашего референс вектора. Если задать vcMinFracLengthF равным 1 то все вектора на вашей картинке станут длинны равной референс вектору (который изображён в правом нижнем углу).
Результат:
Картинка читается и выглядит неплохо, но не хватает похожести на настоящие течения . Вода же не течёт по прямым, поэтому давайте искривим наши векторы. Добавляем в скрипт:
## Change over to curly vectors.
##
vcres.vcMinDistanceF = 0.010
vcres.vcGlyphStyle = "CurlyVector"
vcres.tiMainString = "vcGlyphStyle = 'CurlyVector'"
plot = Ngl.vector_map(wks,u,v,vcres)
Мы немножко меняем минимальное расстояние между векторами - vcres.vcMinDistanceF, чтобы отобразить их больше было на картинке. vcres.vcGlyphStyle , как не трудно догадаться задаёт стиль отображения векторов.
Получаем:
Больше похоже на течения, правда? )
Дальнейшие эксперименты можно проводить долго и я надеюсь что хоть кто ни будь из прочитавших эту статью поделится своими результатами здесь. Также я по мере сил попробую ответить на возникшие вопросы, которые можно оставлять в коментариях.
Под и под конец полностью скрипт:
# Detiled discription in Russian you can find at http://koldunov.net/?p=148
# Script is based on http://www.pyngl.ucar.edu/Examples/Scripts/vector_pop.py
import Ngl
import Nio
import numpy
#
# Open a workstation.
#
rlist = Ngl.Resources()
rlist.wkColorMap = ["White","Black","Tan1","SkyBlue","Red"]
wks_type = "ps"
wks = Ngl.open_wks(wks_type,"vector",rlist)
#
# Open the netCDF file.
#
file = Nio.open_file("./S09_velos_09.nc")
#
# Get the u/v and lat/lon variables.
#
urot = file.variables["SICUO"]
vrot = file.variables["SICVE"]
lat2d = file.variables["lat"]
lon2d = file.variables["lon"]
u = urot[0,0,:,:]
v = vrot[0,0,:,:]
lon = lon2d[:]
lat = lat2d[:]
##
## Now set up a vector resource variable.
##
vcres = Ngl.Resources()
##
## Set coordinate arrays for data. This is necessary in order to overlay
## the vectors on a map.
##
vcres.vfXArray = lon
vcres.vfYArray = lat
vcres.mpProjection = "Orthographic"
vcres.mpFillOn = True
vcres.mpLandFillColor = "Tan1"
vcres.mpOceanFillColor = "SkyBlue"
vcres.mpInlandWaterFillColor = "SkyBlue"
vcres.mpLimitMode = "LatLon"
vcres.mpCenterLonF = -80.
vcres.mpCenterLatF = 55
vcres.mpMinLatF = 60
vcres.mpDataBaseVersion = "MediumRes"
vcres.vcRefAnnoOrthogonalPosF = -0.18 # Move ref anno up into plot.
vcres.tiMainString = "Currents"
plot = Ngl.vector_map(wks,u,v,vcres)
##
## Now we're going to thin the vectors by drawing only every third one.
##
vcres.vfXArray = lon[::3,::3]
vcres.vfYArray = lat[::3,::3]
vcres.tiMainString = "Striding the grid [::3,::3]"
plot = Ngl.vector_map(wks,u[::3,::3],v[::3,::3],vcres)
##
## This time thin the vectors by setting the vcMinDistanceF resource.
## We change the leight of the vectors as well
##
vcres.vfXArray = lon
vcres.vfYArray = lat
vcres.vcRefMagnitudeF = 0.3
vcres.vcMinDistanceF = 0.015
vcres.vcRefLengthF = 0.08
vcres.tiMainString = "vcMinDistanceF = 0.015"
plot = Ngl.vector_map(wks,u,v,vcres)
##
## Make the vectors a little longer.
##
vcres.vcMinFracLengthF = 0.05
vcres.tiMainString = "vcMinFracLenF = 0.05"
plot = Ngl.vector_map(wks,u,v,vcres)
##
## Change over to curly vectors.
##
vcres.vcMinDistanceF = 0.010
vcres.vcGlyphStyle = "CurlyVector"
vcres.tiMainString = "vcGlyphStyle = 'CurlyVector'"
plot = Ngl.vector_map(wks,u,v,vcres)
Классная вещь! SciPy-ем пару раз пользовался, но о такой проге не знал. Надо будет присмотреться к ней поближе.
2fatune
Расскажите о результатах когда присмотритесь )
Хорошая статья! Обязательно поставлю на неё ссылку.
Я вот геоданными не занимаюсь, поэтому ничего толком добавить не могу, но интересно было бы узнать твоё мнение о Basemap в Pylab:
Тоже Python, схожая функциональность (вроде бы тоже умеет визуализировать скалярные и векторные данные в разных геопроекциях). Может пробовал?
2jetxee
Спасибо на добром слове )
Я Basemap не пробовал, но они хотят повторить функционал матлабовского Mapping Toolbox, а на мой вкус он не является идеалом в смысле визуализации данных на карте. Но это не более чем дело вкуса ) По поводу удобства использования и его работы ничего сказать не могу, к сожалению.
Пингль хоть и сделан в климатологической организации, но довольно неплохо справляется с отображением простых графиков, например вот этого , так что думаю и ваши данные ему по зубам )
При этом, конечно, скрипт для отображения данных не маленький. То есть Пингль не годится для быстрого просмотра данных ( типа матлабовского plot(x,y)). Думаю что тот же Basemap (если он действительно по функционалу близок матлабу) больше годится для быстрой, грубой отрисовки. Но создание окончательных, красивых картинок лучше поручать пинглю или GMT.
Хорошо написал, интересно.
Хинт, если в начале подгружать модули таким вот образом:
$: from Ngl import *
$: from math import *
тогда все функции можно будет вызывать пропуская название модуля:
$: wks = open_wks(wks_type,”vector”,rlist)
$: a = sin(b) # вместо a = math.sin(b)
Если кто-то привык к gnuplot, то есть пакет связка python-gnuplot (есть в реаозиториях убунты), с простым синтаксисом.
PS. А с matplotlib есть проблема, кириллица не отображается.
PPS. sudo apt-get install python-*numpy*
Надо всё-таки поправить на ’sudo apt-get install python-numpy’
Понятно, что опечатка. Но некоторые могут просто попробовать скопировать эту строчку в терминал.
И заодно ещё вопрос задам
сейчас наша модель мирового океана использует формат ctl\bin и OpenGrADS для визуализации результатов. Хочу попробовать что-то более красивое. Вот сейчас думаю на счёт PyNGL. А какие ещё “рисовалки” можете порекомендовать/вы используете?
Заранее спасибо.
ps что-то давно ничего не писали в блог)
2dpthebest
Поменял, спасибо, там было аж две опечатки )
Я в общем только PyNGL и использую. Рисовалок достаточно много, начиная от Матлаба, заканчивая всякими самописными велосипедами. В принципе все чем я когда либо пользовался для визуализации описано в этом блоге. В нашей группе в основном используют Ferret, макс-планк для онлайн мониторинга своей модели - cdo, NCAR - естественно NCL. Выбор рисовалки во многом зависит от того сколько времени вы хотите потратить на обучение и какие конкретно задачи нужно решать.
Не писал потому что защита была, а теперь лето - жарко и лениво )