Comparação NCL e GrADS

Neste "post" quero mostrar uma comparação simples entre dois softwares de visualização e processamento de dados amplamente usados na Meteorologia: o NCAR Command Language (NCL) e o Grid Analysis and Display System (GrADS). Ambos são gratuitos e disponíveis na internet nos endereços http://grads.iges.org/grads/grads.html e http://www.ncl.ucar.edu/, respectivamente.

A escolha entre um e outro é definida, principalmente, pelo gosto pessoal ou pela facilidade em se fazer uma determinada tarefa. Geralmente, o motivo preponderante para a definição do pacote a ser usado é a última.

Durante, praticamente, toda a minha vida acadêmica usei o GrADS, basicamente, porque era o único pacote usado na faculdade e na pós-graduação. Os gráficos gerados para o artigo publicado em 2007 na Monthly Weather Review (http://ams.allenpress.com/perlserv/?request=get-abstract&doi=10.1175%2FMWR3302.1) foram gerados com ele.

Recentemente, descobri, por meio do meu orientador de Iniciação Científica, na UFPel, o NCL. Fiquei impressionado com a qualidade dos gráficos e do ambiente de programação proporcionado por ele. Curioso? Faça um tour pelos exemplos do NCL em http://www.ncl.ucar.edu/Applications/.

Voltando ao assunto deste "post", cuja idéia surgiu quando estava fazendo comparações entre os dois, vou mostrar como se calcula a vorticidade relativa e a sua advecção com ambos pacotes e os seus resultados gráficos.

Para começar, veja abaixo o script em GrADS para este propósito:

* abrindo arquivos
'sdfopen uwnd500hPa.nc'
'sdfopen vwnd500hPa.nc'

* definindo janela de plotagem
'set lat -60 -10'
'set lon 280 330'

* vorticidade relativa
'define zeta = hcurl(uwnd.1,vwnd.2)'

* gradiente meridional de zeta
'define dzetax = cdiff(zeta,x)'
'define dy = cdiff(lat,y)*3.1416/180'

* gradiente zonal de zeta
'define dzetay = cdiff(zeta,y)'
'define dx = cdiff(lon,x)*3.1416/180'

* adveccao de zeta
'define adv = -1*( (uwnd.1*dzetax)/(cos(lat*3.1416/180)*dx) + vwnd.2*dzetay/dy )/6.37e6'

* graficos
* PAINEL 1
'set vpage 0 5.5 4.25 8.5'
'set grads off'
'set gxout shaded'
'display zeta'
'set gxout contour'
'display dzetay/(6.37e6*dy)'
'draw title zeta e dzeta/dy'

* PAINEL 2
'set vpage 5.5 11 4.25 8.5'
'set grads off'
'set gxout shaded'
'display zeta'
'set gxout contour'
'display dzetay/(6.37e6*cos(lat*3.1416/180)*dx)'
'draw title zeta e dzeta/dx'

* PAINEL 3
'set vpage 0 5.5 0 4.25'
'set grads off'
'set gxout shaded'
'display zeta'
'set gxout contour'
'display zeta'
'set gxout vector'
'display uwnd.1;vwnd.2'
'draw title zeta e escoamento'

* PAINEL 4
'set vpage 5.5 11 0 4.25'
'set grads off'
'set gxout shaded'
'display zeta'
'set gxout contour'
'display adv'
'draw title zeta e adv(zeta)'

'printim testeGRADS.png white x1000 y1000'

'quit'


Abaixo está o gráfico gerado por este script:



Agora, vamos ao script do NCL para o mesmo cálculo:

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin

; conversão graus -> radianos
g2r = acos(-1.)/180

; abrindo arquivos e carregando variáveis
fu = addfile("uwnd500hPa.nc","r")
fv = addfile("vwnd500hPa.nc","r")
u = short2flt( fu->uwnd(0,0,:,:) )
v = short2flt( fv->vwnd(0,0,:,:) )
lat = fu->lat
dlat = lat(2)-lat(1)
lon = fu->lon
dlon = lon(2)-lon(1)


;------------------ vorticidade relativa
zeta = uv2vr_cfd( u, v, lat, lon, 2 )
copy_VarCoords( u, zeta )

;------------------ gradiente meridional de vorticidade
dzeta_dy = center_finite_diff( zeta(lon|:,lat|:), dlat*g2r*6.37e6, False, 0 )
copy_VarCoords( zeta(lon|:,lat|:), dzeta_dy)

;------------------ gradiente zonal de vorticidade
dzeta_dx = new( dimsizes(u), typeof(u), u@_FillValue )
do k = 0, dimsizes( lat )-1
dx = cos( lat(k)*g2r )*dlon*g2r*6.37e6
dzeta_dx(k,:) = center_finite_diff( zeta(k,:), dx, False, 0 )
end do
copy_VarCoords( zeta, dzeta_dx )

;------------------ adveccao de vorticidade relativa
adv_zeta = -1*( u(lat|:,lon|:)*dzeta_dx(lat|:,lon|:) + \
v(lat|:,lon|:)*dzeta_dy(lat|:,lon|:) ) * 1e10
copy_VarCoords( zeta, adv_zeta )

; janela de plotagem
minlat = -60
maxlat = -10
minlon = 280
maxlon = 330

;------------------ graficos
wks = gsn_open_wks( "ps", "advVorticidadeNCL" )
gsn_define_colormap(wks,"BlWhRe")

; recursos para os gráficos
res = True
res@gsnDraw = False
res@gsnFrame = False
res@gsnMaximize = True
res@gsnAddCyclic = False
res@mpMinLonF = minlon
res@mpMaxLonF = maxlon
res@mpMinLatF = minlat
res@mpMaxLatF = maxlat
res@cnFillOn = True
res@lbLabelBarOn = False
res@gsnSpreadColors = True

res2 = True
res2@gsnAddCyclic = False

; arranjo para gráficos
plot = new( 4, graphic )

; PAINEL 1
res@gsnLeftString = "zeta e dzeta/dy"
plot(0) = gsn_csm_contour_map_overlay( wks, \
zeta({lat|minlat:maxlat},{lon|minlon:maxlon}), \
dzeta_dy({lat|minlat:maxlat},{lon|minlon:maxlon}), \
res, res2 )

; PAINEL 2
res@gsnLeftString = "zeta e dzeta/dx"
plot(1) = gsn_csm_contour_map_overlay( wks, \
zeta({lat|minlat:maxlat},{lon|minlon:maxlon}), \
dzeta_dx({lat|minlat:maxlat},{lon|minlon:maxlon}), \
res, res2 )

; PAINEL 3
res@gsnLeftString = "zeta e escoamento"
res@gsnScalarContour = True ; para plotar contornos e vetores separadamente
plot(2) = gsn_csm_vector_scalar_map( wks, \
u({lat|minlat:maxlat},{lon|minlon:maxlon}), \
v({lat|minlat:maxlat},{lon|minlon:maxlon}), \
zeta({lat|minlat:maxlat},{lon|minlon:maxlon}), \
res )

; PAINEL 4
res@gsnLeftString = "zeta e adv(zeta)"
plot(3) = gsn_csm_contour_map_overlay( wks,
zeta({lat|minlat:maxlat},{lon|minlon:maxlon}), \
adv_zeta({lat|minlat:maxlat},{lon|minlon:maxlon}), \
res, res2 )

; plotando os quatro paineis numa pagina
gsn_panel( wks, plot, (/2,2/), False )

; conversão PS->PNG
system("convert -density 300 -geometry 1000x1000 -trim advVorticidadeNCL.ps advVorticidadeNCL.png")

end


Abaixo está o gráfico resultante:



Os gráficos gerados são aqueles padrões para ambos pacotes. Tentei interfirir o mínimo possível nos mesmo, para que a comparação fosse a mais justa possível.

Os recursos gráficos para os campos plotados com o NCL (definidos por meio dos atributos passados as variáveis res e res2; uso do símbolo "@") são aqueles estritamente necessários para a definição do domínio de plotagem (mp*LonF e mp*LatF), plotagem de campos sobrepostos (gsnFrame e gsnDraw), plotagem de dados não globais (gsnAddCyclic), uso de cores (cnFillOn e gsnSpreadColors) e para plotagem de vetores e contornos (gsnScalarContour). A única linha que muda o gráfico padrão é a

gsn_define_colormap(wks,"BlWhRe")

que escolhe um outro conjunto de cores para campos com preenchimento de cores.

A qualidade superior dos gráficos gerados pelo NCL é facilmente percebida, mas é na parte de programação que vejo a grande vantagem oferecida pelo NCL em comparação com o GrADS. A possibilidade de se lidar com vetores, matrizes e arranjos de variadas dimensões é um forte exemplo desta vantagem.

Ambos pacotes têm funções para o cálculo da vorticidade e para o cálculo de diferenças finitas.

O único ponto no qual o NCL deixa a desejar (ao meu ver) é a perda das informações de coordenadas (os metadados) quando se faz qualquer operação matemática com uma variável, o que não ocorre no GrADS. Suponha que tenho o campo da componente meridional do vento, dada por v[lat, lon], dada em m/s. No GrADS, a operação

vKmH = v*3.6

pode ser plotada diretamente. Já no NCL, preciso copiar as informações das coordenadas espaciais para a nova variável antes de plotá-la

copy_VarCoords( v, vKmH )

Este procedimento foi usado várias vezes no script NCL acima.

O GrADS também tem a facilidade de gerar gráficos em GrADS Metafile (MF), PNG, JPEG e GIF diretamente, enquanto o NCL gera seus gráficos, diretamente, em PS, EPS, EPSI, PDF e NCGM. O GrADS oferece um programa externo (gxps) para a conversão de MF para PS.

A figura PNG foi gerada com o GrADS com a função printim. O arquivo PS do NCL foi convertido para PNG com o programa externo convert, disponibilizado pelo pacote ImageMagick, no linux.

Este é um pequeno exemplo de comparação entre os dois pacotes. Cada um tem as suas vantagens e desvantagens. Não há "o melhor" pacote, mas sim aquele mais adequado às preferências e necessidades pessoais. Faça a sua escolha, compare! Se você usou apenas um deles até hoje, teste o outro!

Em "posts" futuros colocarei mais exemplos de comparações entre os dois pacotes.

Abraços,

INFO: Os dados usados para a geração destes gráficos são da Reanálise 2 do NCEP (http://www.cdc.noaa.gov/data/gridded/data.ncep.reanalysis2.pressure.html). São dados médios do dia 01/01/1979 para as duas componentes do vento: u e v, em 500 hPa.