Одновременное отображение скалярных и векторных данных в PyNGL

Задача: построить высоту поверхности моря и наложить на неё геострофические течения.
?нструменты: Python, PyNGL

В предыдущих постах мы рисовали векторные и скалярные данные. Теперь попробуем их объединить.

Загружать данные на этот раз будем из открытых, постоянно обновляемых источников, так что каждый сможет ежедневно следить за тем как там поживает Гольфстрим, какие вихри от него отделились и куда эти вихри движутся. Выглядеть это в итоге будет примерно вот так:

pyngl_small1.png

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

Сначала приведу скрипт целиком:

import Ngl
import Nio
import numpy
import sys


try:
    h_file = sys.argv[1];  uv_file = sys.argv[2]
except:
    print "Usage:",sys.argv[0], "h_file uv_file"
    sys.exit(1)
   

file_h = Nio.open_file(h_file)
file_uv = Nio.open_file(uv_file)


height_nc = file_h.variables["Grid_0001"]
lat_nc  = file_h.variables["NbLatitudes"]
lon_nc = file_h.variables["NbLongitudes"]
u = file_uv.variables["Grid_0001"]
v = file_uv.variables["Grid_0002"]


lat = lat_nc[:]
lon = lon_nc[:]
height = height_nc[:]
uu = u[:]
vv = v[:]


latt = numpy.zeros((1080,915))
lonn = numpy.zeros((1080,915))

for i in range(1080):
    latt[i,:] = lat[:]
   
for i in range(915):
    lonn[:,i] = lon[:]
   
   
rlist            = Ngl.Resources()
rlist.wkColorMap = 'posneg_1'
wks_type = "ps"
wks = Ngl.open_wks(wks_type,"ssh",rlist)

ws_id = Ngl.get_workspace_id()
rlist = Ngl.Resources()
rlist.wsMaximumSize = 33554432
Ngl.set_values(ws_id,rlist)

##
## Set up resource lists for map, contour, and vector
## plots.
##
mapres = Ngl.Resources()
cnres  = Ngl.Resources()
vecres = Ngl.Resources()

##
##  Set the scalarfield missing value if height_nc has one specified.
##
if hasattr(height_nc,"_FillValue"):
    cnres.sfMissingValueV = float(height_nc._FillValue[0])
    vecres.vfMissingUValueV = float(height_nc._FillValue[0])
    vecres.vfMissingVValueV = float(height_nc._FillValue[0])

##
## Turn off nglDraw and nglFrame because we don't want to draw all
## these plots until they are all overlaid on the map plot.
##

mapres.nglDraw  = False
mapres.nglFrame = False
cnres.nglDraw  = False
cnres.nglFrame = False
vecres.nglDraw  = False
vecres.nglFrame = False



## Map related resources
mapres.mpProjection          = "Mercator"
mapres.mpDataBaseVersion     = "MediumRes"
mapres.mpLimitMode           = "LatLon"
mapres.mpMinLonF             = -90
mapres.mpMaxLonF             = -45
mapres.mpMinLatF             = 22
mapres.mpMaxLatF             = 45
mapres.mpFillOn              = True  

igray = Ngl.new_color(wks,0.7,0.7,0.7)
mapres.mpFillColors = [0,-1,igray,-1]

mapres.mpLandFillColor       = igray
mapres.mpOceanFillColor      = -1
mapres.mpGridAndLimbOn       = False


## Contours releted resources

cnres.sfXArray        = lonn[:]
cnres.sfYArray        = latt[:]

cnres.cnFillOn             = True
cnres.cnFillDrawOrder       = "Predraw"
cnres.cnLineLabelsOn        = False
cnres.cnLinesOn                 = False
cnres.nglSpreadColorStart  = 5
cnres.nglSpreadColorEnd      = -2
cnres.cnLineDrawOrder      = "Predraw"
cnres.cnLevelSelectionMode = "ExplicitLevels" # Define own levels.
cnres.cnLevels             = numpy.arange(45.,215.,20)

## Legend
cnres.lbTitleString  = height_nc.__dict__['long_name']+', '+height_nc.__dict__['units']
cnres.lbTitleFontHeightF        = 0.018
cnres.lbLabelFontHeightF        = 0.012
cnres.lbTitleOffsetF            = -0.27
cnres.lbBoxMinorExtentF         = 0.15
cnres.pmLabelBarOrthogonalPosF  = -0.01
cnres.lbOrientation             = "Horizontal"

## Vector related resources

vecres.vfXArray = lonn[:]
vecres.vfYArray = latt[:]

#vecres.vcRefAnnoOrthogonalPosF = -0.22
vecres.vcRefAnnoOn    = False
vecres.vcRefMagnitudeF    = 30.0
vecres.vcRefLengthF       = 0.08
vecres.vcMinDistanceF     = 0.01
vecres.vcGlyphStyle       = "CurlyVector"
vecres.vcVectorDrawOrder = "Predraw"

#
# Create the various plots. They will not get drawn because
# nglDraw is set to False for all of them.
#

map_plot          = Ngl.map(wks,mapres)
vector_plot       = Ngl.vector(wks,uu,vv,vecres)
line_contour_plot = Ngl.contour(wks,height[:,:],cnres)

#
# Overlay everything on the map plot.
#

Ngl.overlay(map_plot,line_contour_plot)
Ngl.overlay(map_plot,vector_plot)

#
# Change the title.
#
srlist = Ngl.Resources()
srlist.tiMainString = "Data by "+file_h.__dict__['CreatedBy']+". "+file_h.__dict__['CreatedOn']
#srlist.tiMainOffsetYF    = 0.025
srlist.tiMainFontHeightF = 0.024
Ngl.set_values(map_plot,srlist)

#
# Draw the map plot, which now contains the vectors and
# filled/line contours.
#

Ngl.maximize_plot(wks,map_plot)    # Maximize size of plot in frame.
Ngl.draw(map_plot)
Ngl.frame(wks)

Много картинок по ходу не будет, просто объясню некоторые части.

Данные мы будем брать с сайта AVISO, а конкретно нам нужны Maps of Absolute Dynamic Topography & absolute geostrophic velocities. Первые берём отсюда, вторые отсюда. Сервера у них немного тормозные. Файлы по полтора и три метра соответственно, расширение *.nc.gz. Поскольку денег мы не платим, данные которые там выложены с задержкой 30 дней, но для научных целей этого вполне достаточно.

?звлекаем данные из архива и передаём имена файлов в качестве аргументов скрипту, например так:

python ssh.py madt_oer_merged_h_21247.nc madt_oer_merged_uv_21247.nc

Общий принцип работы скрипта следующий.

  • создаём отдельно переменные ресурсов для различных слоёв — карты, изолиний и векторов
  • задаём нужные нам ресурсы (то есть определяем их свойства)
  • виртуально строим эти карты (не выводя ничего)
  • сводим три слоя вместе и выводим это дело в файл

Теперь посмотрим на некоторые части скрипта.
Сначала немного преобразуем входные данные:

latt = numpy.zeros((1080,915))
lonn = numpy.zeros((1080,915))

for i in range(1080):
    latt[i,:] = lat[:]
   
for i in range(915):
    lonn[:,i] = lon[:]

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

и вектор долгот lon
lon.png

?х них нам надо сделать две одинаковых по размерам матрицы:
latt
llat.png

и lonn
llon.png

Конечно наши матрицы немножко больше и вследствие того что расчёты ведутся на проекции Меркатора расстояние между точками сетки по широте различное (по долготе оно одинаковое и равняется примерно 1/3 градуса).

Далее несколько загадочных строчек:

ws_id = Ngl.get_workspace_id()
rlist = Ngl.Resources()
rlist.wsMaximumSize = 33554432
Ngl.set_values(ws_id,rlist)

Этот код советуют вставлять в случае затыка который произошёл со мной:

fatal:ContourPlotDraw: Workspace reallocation would exceed maximum size 16777216

Памяти рабочей станции не хватает, и лечится всё таким образом. В подробности я не вдавался )

Следующие строчки жизненно важны:

##
##  Set the scalarfield missing value if height_nc has one specified.
##
if hasattr(height_nc,"_FillValue"):
    cnres.sfMissingValueV = float(height_nc._FillValue[0])
    vecres.vfMissingUValueV = float(height_nc._FillValue[0])
    vecres.vfMissingVValueV = float(height_nc._FillValue[0])

Они сообщают что мы будем считать отсутствующими значениями. Значения эти есть в .nc файле, и когда мы ассоциируем переменную height_nc с переменной Grid_0001, то можем добраться до её так называемых Missing values значения при помощи конструкции height_nc._FillValue. Для u и v можно было тянуть значения из соответствующей переменной, но поскольку они одинаковые, заморачиваться я не стал 🙂 . Да, и if тут в принципе не обязателен.

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

Создаём несколько слоёв:

map_plot          = Ngl.map(wks,mapres)
vector_plot       = Ngl.vector(wks,uu,vv,vecres)
line_contour_plot = Ngl.contour(wks,height[:,:],cnres)

Но они нарисованы не будут, а будут помещены в переменные.

Теперь последовательно накладываем на слой карты слои изолиний и векторов:

#
# Overlay everything on the map plot.
#

Ngl.overlay(map_plot,line_contour_plot)
Ngl.overlay(map_plot,vector_plot)

Почему то наложить map_plot на line_contour_plot, то есть что то типа Ngl.overlay(line_contour_plot,map_plot) у меня не получилось. Пытался я делать это потому что сначала часть векторов наезжала на сушу. Но оказалось что это вполне адеквано контролируется ресурсами vcVectorDrawOrder, cnLineDrawOrder и cnFillDrawOrder.

Добавляем заголовок:

srlist = Ngl.Resources()
srlist.tiMainString = "Data by "+file_h.__dict__['CreatedBy']+". "+file_h.__dict__['CreatedOn']
#srlist.tiMainOffsetYF    = 0.025
srlist.tiMainFontHeightF = 0.024
Ngl.set_values(map_plot,srlist)

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

Ну и под конец рисуем всё это дело в выходной файл.
Думаю наладить автоматическое производство этих карт и постить их либо тут, либо на oceanographers 🙂

PyNGL H and currents