domingo, 28 de fevereiro de 2010

Funções externas no NCL

O NCL possui uma grande variedade de funções e subrotinas, como pode ser visto em http://www.ncl.ucar.edu/Document/Functions/, que mostra as funções do NCL organizadas por categorias. Mesmo com todas estas funções alguma operação bem específica pode não possuir uma função ou subrotina no NCL e você deverá produzir um programa para que ela seja realizada.

Dependendo da complexidade da operação, a utilização das funções da linguagem de script do NCL para criar este programa pode torná-lo ineficaz, levando muito tempo para completar a operação ou mesmo não conseguindo realizá-la. Nestas situações é melhor utilizar uma subrotina ou função externa, criada em outra linguagem de programação.

Um exemplo de uma operação não existente no NCL é a de adicionar linhas a um arquivo ASCII já existente. Aqui vou mostrar como é possível adicionar uma função deste tipo ao NCL por meio da linguagem Fortran 90. O NCL possibilita a utilização de subrotinas ou funções escritas em Fortran (77 ou 90) e em C nos seus scripts. Maiores informações sobre este procedimento podem ser obtidas em http://www.ncl.ucar.edu/Document/Manuals/Ref_Manual/NclExtend.shtml.

Apresento abaixo os três passos necessários para a utilização da subrotina escrita em Fortran 90.

Passo 1) Escrever a subrotina em Fortran 90, que será salva num arquivo chamado adicionaLinha.f90:

subroutine addline(arquivo, dimLinha, fmt1, fmt2, linha)

implicit none


integer,parameter :: unidES=30
logical :: existeArq
integer,intent(in) :: dimLinha,fmt1,fmt2
character(len=255),intent(in) :: arquivo
real,dimension(dimLinha),intent(in) :: linha

inquire(file=trim(adjustl(arquivo)),exist=existeArq)

if (.NOT.existeArq) then

print *,"Aviso: arquivo nao existe, arquivo sendo criado!"
open(unidES,file=trim(adjustl(arquivo)),status='new')
write(unidES,'(F.)') linha(1:dimLinha)

else

open(unidES,file=trim(adjustl(arquivo)),access='append')
write(unidES,'(F.)') linha(1:dimLinha)

end if

close(unidES)

end subroutine addLine


Passo 2) Criar um arquivo auxiliar que permite o NCL entender a sintaxe do programa em Fortran 90 apresentado no Passo 1. Este arquivo contém uma declaração da subrotina em Fortran 90 na sintaxe do Fortran 77 (Note a diferença entre as declarações das variáveis da subrotina.). Neste exemplo, o arquivo auxiliar se chamará adicionaLinha.stub:

C NCLFORTSTART
subroutine addline(arquivo,dimLinha,fmt1,fmt2,linha)
character*(*) arquivo
integer dimLinha,fmt1,fmt2
real linha(dimLinha)
C NCLEND

Passo 3) Criar um arquivo objeto compartilhado (.so) usando o comando WRAPIT (http://www.ncl.ucar.edu/Document/Tools/WRAPIT.shtml), que vem com o próprio NCL. O comando WRAPIT possui várias opções. Neste exemplo, uso uma em especial, a -in, que manda o WRAPIT usar o compilador Fortran da Intel (http://software.intel.com/en-us/articles/non-commercial-software-download/) para criar o objeto compartilhado (Por favor, veja nas notas incluídas ao final deste texto do porquê da utilização deste compilador), sendo executado da seguinte forma:

WRAPIT -in adicionaLinha.stub adicionaLinha.f90

Ao final da execução desta linha de comando um arquivo adicionaLinha.so será criado. É este arquivo que pode ser usado em scripts do NCL para executar a operação programada na subrotina escrita em Fortran 90. Abaixo passo um script em NCL para testar esta função externa:

external ADDLINE "./adicionaLinha.so"

begin

arq = "teste.txt"
a = (/1.,2.,3.,4.,5.,6.,7.,8.,9.,10./)
dimA = dimsizes(a)
fmt1 = 5
fmt2 = 1

ADDLINE::addline(arq,dimA,fmt1,fmt2,a)


end

Pronto! Uma nova função foi adicionada ao NCL graças a sua interação com outras linguagens de programação, como o Fortran 90. A execução do script em NCL escrito acima resultará num arquivo chamado teste.txt. Caso este arquivo já exista no diretório de execução do script será adicionada uma nova linha ao final do arquivo. Faça um teste, execute mais de uma vez este mesmo script!

É isso aí, espero que esta dica seja útil.



NOTA) Por que é usado o compilador Fortran da Intel?

O compilador Fortran da Intel (de uso não comercial) é usado neste exemplo por possuir algumas características do Fortran 95, como a possibilidade de criação de Expressões Variáveis de Formato (http://www.intel.com/software/products/compilers/docs/flin/main_for/lref_for/source_files/pghvarfm.htm). Isto é aplicado na subrotina escrita em Fortran 90, nas linhas com o comando write:

write(unidES,'(F.)') linha(1:dimLinha)

Nesta linha, as variáveis dimLinha, fmt1 e fmt2 são usadas para controlar o formato dos valores da variável numérica linha que serão impressos no arquivo ASCII. Essa característica proporciona a criação de uma subrotina mais flexível.

terça-feira, 13 de outubro de 2009

Acentos no NCL

Não há no NCL nenhuma fonte com caracteres acentuados e nem a interpretação automática de strings com tais caracteres. Isso atrapalha os usuários cujas línguas nativas exigem estes caracteres. É possível escrever palavras com acentos no NCL, mas é um pouco chato.

Os caracteres acentuados são gerados por meio de Códigos de Função, que permite controlar vários aspectos de um texto. Veja a descrição completa destes códigos aqui.

Para facilitar a minha vida, criei um arquivo que armazena separadamente os códigos para os caracteres acentuados. Os nomes das variáveis do tipo string, que contêm os códigos, seguem os nomes da tabela de acentos do HTML. Veja abaixo os caracteres acentuados e algumas palavras em português, mostrando como o NCL os plota.


Abaixo está o script que gerou esta saída:


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

begin

; tabela de acentos
Agrave = "A~H-15V6F35~A~FV-6H3~" ; À
agrave = "a~H-13V2F35~A~FV-2H3~" ; à
Aacute = "A~H-15V6F35~B~FV-6H3~" ; Á
aacute = "a~H-13V2F35~B~FV-2H3~" ; á
Acirc = "A~H-15V6F35~C~FV-6H3~" ; Â
acirc = "a~H-13V2F35~C~FV-2H3~" ; â
Atilde = "A~H-15V6F35~D~FV-6H3~" ; Ã
atilde = "a~H-13V2F35~D~FV-2H3~" ; ã
Auml = "A~H-15V6F35~H~FV-6H3~" ; Ä
auml = "a~H-13V2F35~H~FV-2H3~" ; ä

Egrave = "E~H-15V6F35~A~FV-6H3~" ; È

egrave = "e~H-13V2F35~A~FV-2H3~" ; è
Eacute = "E~H-15V6F35~B~FV-6H3~" ; É
eacute = "e~H-13V2F35~B~FV-2H3~" ; é
Ecirc = "E~H-15V6F35~C~FV-6H3~" ; Ê
ecirc = "e~H-13V2F35~C~FV-2H3~" ; ê
Euml = "E~H-15V6F35~H~FV-6H3~" ; Ë
euml = "e~H-13V2F35~H~FV-2H3~" ; ë

Igrave = "I~H-10V6F35~A~FV-6H3~" ; Ì

igrave = "i~H-10V2F35~A~FV-2H3~" ; ì
Iacute = "I~H-08V6F35~B~FV-6H3~" ; Í
iacute = "i~H-08V2F35~B~FV-2~" ; í
Icirc = "I~H-09V6F35~C~FV-6H3~" ; Î
icirc = "i~H-09V2F35~C~FV-2H3~" ; î
Iuml = "I~H-09V6F35~H~FV-6H3~" ; Ï
iuml = "i~H-09V2F35~H~FV-2H3~" ; ï

Ograve = "O~H-15V6F35~A~FV-6H3~" ; Ò

ograve = "o~H-13V2F35~A~FV-2H3~" ; ò
Oacute = "O~H-15V6F35~B~FV-6H3~" ; Ó
oacute = "o~H-13V2F35~B~FV-2H3~" ; ó
Ocirc = "O~H-16V6F35~C~FV-6H3~" ; Ô
ocirc = "o~H-14V2F35~C~FV-2H3~" ; ô
Otilde = "O~H-15V6F35~D~FV-6H3~" ; Õ
otilde = "o~H-13V2F35~D~FV-2H3~" ; õ
Ouml = "O~H-16V6F35~H~FV-6H3~" ; Ä
ouml = "o~H-14V2F35~H~FV-2H3~" ; ä

Ugrave = "U~H-15V6F35~A~FV-6H3~" ; Ù
ugrave = "u~H-13V2F35~A~FV-2H3~" ; ù
Uacute = "U~H-13V6F35~B~FV-6H3~" ; Ú
uacute = "u~H-13V2F35~B~FV-2H3~" ; ú
Ucirc = "U~H-15V6F35~C~FV-6H3~" ; Û
ucirc = "u~H-13V2F35~C~FV-2H3~" ; û
Uuml = "U~H-15V6F35~H~FV-6H3~" ; Ü
uuml = "u~H-13V2F35~H~FV-2H3~" ; ü

Cedil = "C~H-15F35~K~FH2~" ; Ç
cedil = "c~H-13F35~K~FH2~" ; ç

Ntilde = "N~H-15V6F35~D~FV-6H3~" ; Ñ
ntilde = "n~H-13V2F35~D~FV-2H3~" ; ñ

; ambiente grafico
wks = gsn_open_wks("ps","acentos")

; recursos do texto

txres = True
txres@txFontHeightF = 0.03

; plotando caracteres acentuados
txres@txJust = "CenterCenter"
titulo1 = "Acentua"+cedil+atilde+"o com o NCL"
gsn_text_ndc(wks,titulo1,.5,.95,txres)

txres@txJust = "CenterLeft"
titulo2 = "Caracteres acentuados:"
gsn_text_ndc(wks,titulo2,0.,.85,txres)

texto = Agrave+" "+agrave+" "+Aacute+" "+\
aacute+" "+Acirc+" "+acirc+" "+\
Atilde+" "+atilde+" "+ \
Auml+" "+auml
gsn_text_ndc(wks,texto,.2,.8,txres)

texto = Egrave+" "+egrave+" "+Eacute+" "+\
eacute+" "+Ecirc+" "+ecirc+" "+\
Euml+" "+euml
gsn_text_ndc(wks,texto,.2,.75,txres)

texto = Igrave+" "+igrave+" "+Iacute+" "+\
iacute+" "+Icirc+" "+icirc+" "+\
Iuml+" "+iuml
gsn_text_ndc(wks,texto,.2,.7,txres)

texto = Ograve+" "+ograve+" "+Oacute+" "+\
oacute+" "+Ocirc+" "+ocirc+" "+\
Otilde+" "+otilde+" "+\
Ouml+" "+ouml
gsn_text_ndc(wks,texto,.2,.65,txres)

texto = Ugrave+" "+ugrave+" "+Uacute+" "+\
uacute+" "+Ucirc+" "+ucirc+" "+\
Uuml+" "+uuml
gsn_text_ndc(wks,texto,.2,.6,txres)

texto = Cedil+" "+cedil+" "+Ntilde+" "+ntilde
gsn_text_ndc(wks,texto,.2,.55,txres)

; algumas palavras acentuadas

titulo3 = "Algumas palavras:"
gsn_text_ndc(wks,titulo3,0.,.45,txres)

txres@txFontHeightF = 0.025
texto = "ver"+atilde+"o - inst"+aacute+\
"vel - hist"+oacute+"rico - mat"+eacute+"ria"
gsn_text_ndc(wks,texto,.1,.4,txres)

texto = "precipita"+cedil+atilde+\
"o - INFORMA"+Cedil+Otilde+"ES - tr"+ecirc+"s"
gsn_text_ndc(wks,texto,.1,.35,txres)

texto = "dire"+cedil+otilde+\
"es - El-Ni"+ntilde+"o - LA-NI"+Ntilde+"A"
gsn_text_ndc(wks,texto,.1,.3,txres)

texto = "got"+iacute+"culas - pol"+ecirc+\
"mica - Amaz"+ocirc+"nia"
gsn_text_ndc(wks,texto,.1,.25,txres)

texto = "mec"+acirc+"nica - portugu"+ecirc+\
"s - ci"+ecirc+"ncia atmosf"+eacute+"rica"
gsn_text_ndc(wks,texto,.1,.2,txres)

texto = Aacute+"reas - din"+acirc+"mica - Sa"+\
uacute+"de"
gsn_text_ndc(wks,texto,.1,.15,txres)

frame(wks)

end


Neste exemplo, os caracteres acentuados estão no script que os usam, mas podem ser colocados num arquivo .ncl separado que pode ser carregado com o comando load.

A fonte de caracteres 35 contém os caracteres de acentuação, sendo usada, extensivamente, na criação dos caracteres acentuados. Vou analisar a construção de um dos caracteres, sendo direta a extensão desta explicação aos demais. Vejamos o caractere "á", definido no script acima como

aacute = "a~H-13V2F35~B~FV-2H3~"

A explicação se dará quase caractere por caractere:

a é o próprio caractere a
~ é o código de função de texto, que pode ser modificado no arquivo .hluresfile.
H-13 move a posição do texto atual 13 pontos na horizontal, à esquerda (negativo).
V2 move a posição do texto atual 2 pontos na vertical, para cima.
F35 usa a fonte de caracteres 35.
B usa o caractere da fonte de caracteres 35.
F volta a usar a fonte de caracteres anterior.
V-2 move a posição do texto atual 2 pontos na vertical, para baixo.
H3 move a posição do texto na horizontal, três pontos à direita.

Em suma, os caracteres acentuados são montados manualmente, por meio de movimentos do cursor de impressão de caractere para a direita, esquerda, para cima e para baixo.

Quem usou uma máquina de escrever já percebeu o esquema de acentuação, ou seja, (i) escreve o caractere, (ii) volta o cursor para a posição do caractere escrito, (iii) eleva ou abaixa o cursor, (iv) escreve o sinal de acentuação e (v) faz o caminho de volta para continuar a escrever o texto.

Note que os códigos de função de texto ficam sempre entre ~, o sinal indicador do uso destes códigos.

Como comentei no início deste post, a criação dos caracteres acentuados é chata, mas é simples, como se pode ver acima. Uma vez feita, é só sair usando ;).

Abraços e espero que seja útil.

quinta-feira, 8 de outubro de 2009

Unindo figuras com o ImageMagick

Imaginemos a seguinte situação: inúmeras figuras individuais foram geradas por um script e/ou programa, mas somente depois do demorado processamento chegamos à conclusão que o ideal era agrupar as figuras de um modo melhor. O que fazer? Processamos tudo de novo, gerando as figuras da maneira que queremos?

A resposta: sim e não. Sim, se temos tempo disponível para alterar o script e processar os dados novamente para gerar as novas figuras. Não, se o tempo está escasso.

Vamos supor que o tempo está escasso e devemos nos virar com as figuras que temos. O passo óbvio a seguir seria usar um editor de imagens ou editor de texto e agrupar as figuras que queremos. Infelizmente, dá trabalho. Cliques com mouse, seleção das figuras etc. Que tal colocarmos as figuras que queremos agrupar numa linha de comando e, como resultado, tê-las unidas da maneira que queremos.

A solução: o ImageMagick, um pacote de comandos para manipulação de imagens. Entre os vários comandos disponíveis neste pacote, há um chamado montage usado para se criar montagens com figuras.

Vamos aos exemplos práticos. Para isso, vamos supor que as figuras geradas com o NCL no post anterior, ao invés de agrupadas como estão, foram geradas individualmente e vamos agrupar as quatro primeiras daquelas figuras, cujos nomes são: leituraBinarioETA-0.png, leituraBinarioETA-1.png, leituraBinarioETA-2.png e leituraBinarioETA-3.png

Exemplo 1: montagem simples e direta das quatro figuras

montage leituraBinarioETA-0.png leituraBinarioETA-1.png leituraBinarioETA-2.png leituraBinarioETA-3.png montagem1.png

o resultado:


Este é o resultado padrão, sem o uso de opções para o comando montage. Clique na figura e veja o real tamanho desta figura.


Exemplo 2: mesma montagem, mas com o tamanho original das figuras e espaços entre elas

montage -geometry +50+50 leituraBinarioETA-0.png leituraBinarioETA-1.png leituraBinarioETA-2.png leituraBinarioETA-3.png montagem2.png

o resultado:


A opção -geometry +50+50 adiciona espaços na horizontal e na vertical entre as figuras. Clique na figura e veja o tamanho real dela.


Exemplo 3: montando uma coluna com as quatro figuras

montage -tile 1x4 -geometry +0+50 leituraBinarioETA-0.png leituraBinarioETA-1.png leituraBinarioETA-2.png leituraBinarioETA-3.png montagem3.png

o resultado:


Aqui, a opção -tile 1x4 diz ao comando montage para criar uma "figura coluna" com as quatro figuras passadas ao comando.


Exemplo 4: montando uma linha com as quatro figuras

montage -tile 4x1 -geometry +50+0 leituraBinarioETA-0.png leituraBinarioETA-1.png leituraBinarioETA-2.png leituraBinarioETA-3.png montagem4.png

o resultado:



Agora, a mesma opção -tile é usada para criar uma "figura linha" com as quatro figuras.


Exemplo 5: igual ao Exemplo 2, mas com a inclusão de texto para cada uma das figuras

montage -geometry +50+50 -label 'a' leituraBinarioETA-0.png -label 'b' leituraBinarioETA-1.png -label 'c' leituraBinarioETA-2.png -label 'd' le ituraBinarioETA-3.png montagem5.png

o resultado:


Neste exemplo, a inclusão de texto logo abaixo das figuras, passadas ao comando pela opção -label.


Mais exemplos e opções para o comando montage pode ser encontrado aqui. Com estes exemplos vimos como é fácil unir figuras de várias formas com o comando montage. Há opções para colocar bordas, usar informações dos próprios arquivos para o texto que acompanha as figuras e muito mais.

Espero que seja útil.

Abraços.

PS. (09/10/2009): quando escrevi este post eu não sabia como controlar o tamanho do texto dos 'labels' atribuídos às figuras. Veja que no Exemplo 5, os 'labels' são extremamente pequenos. Isto é feito com a opção -pointsize tam, sendo tam um número indicando o tamanho do texto. Para o Exemplo 5, um tamanho igual a 50 produz um bom resultado.

segunda-feira, 5 de outubro de 2009

Leitura de arquivos binários com o NCL

O NCL permite a leitura e escrita de arquivos binários, além de outros tipos de arquivos. Há funções para ler arquivos binários gerados tanto por C quanto por Fortran. Como não uso o C, vou mostrar um exemplo com uma das funções para binários gerados por Fortran.

A função usada neste exemplo chama-se "fbindirread", que lê arquivos binários de acesso direto. O script que apresento neste post lê um arquivo binário gerado pelo pós processamento do modelo regional ETA. Baixe o arquivo binário (download) para testar o script NCL apresentado mais abaixo.

O arquivo de controle (CTL) para visualização no GrADS é mostrado abaixo:

DSET ^Indices2006032700+2006032718.bin
UNDEF -9999.
TITLE Indices de Tempestades Severas
XDEF 140 LINEAR -56.0 0.1000
YDEF 85 LINEAR -26.5 0.1000
ZDEF 1 LEVELS 1000
TDEF 1 LINEAR 00Z27Mar2006 1hr
VARS 17
ncl 0 99 Nivel de Cond. Levantamento (m)
nce 0 99 Nivel de Conv. Espontanea (m)
ne 0 99 Nivel de Equilibrio (m)
thetab 0 99 Max dif. Temp. Bulbo umido (K)
li 0 99 Lifted Index (degree)
cape 0 99 CAPE (J/kg)
mcape 0 99 CAPE (J/kg)
cin 0 99 CINE (J/kg)
shr37 0 99 Total 6h Precip. (m)
ustrm 0 99 Desl. Esp. da Tempestade (m/s)
vstrm 0 99 Desl. Esp. da Tempestade (m/s)
heli 0 99 Helicidaade Rel. a Tempestade (m^2/s^2)
nrv 0 99 Numero de Richardson Volumetrico ()
dnrv 0 99 Den. de Num.de R. Volumetrico (m^2/s^2)
dnrv2km 0 99 Den. de Num.de R. V. Cam. 2km (m^2/s^2)
ieh 0 99 Indice de Energia-Helicidade ()
sup 0 99 Parametro de Super Celula ()
ENDVARS

E agora, o script NCL que le o arquivo binário e plota suas variáveis:

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


; ============== LEITURA DOS DADOS BINARIOS ============
arquivo = "Indices2006032700+2006032718.bin"


; dimensoes das variaveis

; informação obtida do CTL
lon = fspan( -56., -56+139*0.1, 140 )
lat = fspan( -26.5, -26.5+84*0.1, 85 )


; lendo todos os 17 campos do arquivo
saidasETA = fbindirread( arquivo, 0, (/17,85,140/), "float" )

; ==== ATRIBUICAO DE COORDENADAS E OUTRAS INFORMACOES =====
; atribuindo coordenadas ao campo lido

saidasETA!0 = "campos"

saidasETA!1 = "lat"

saidasETA&lat = lat
saidasETA&lat@units = "degree_north"
saidasETA!2 = "lon"

saidasETA&lon = lon
saidasETA&lon@units = "degree_east"


; nome dos campos contidos no arquivo binario
nomeCampos = (/"NCL","NCE","NE","Max. Dif. T~B~W~N~","LI","CAPE", \
"MCAPE", "CINE", "Precip. em 6h", "Desloc. U da Tempestade", \
"Desloc. V da Tempestade", "Helic. Rel. Tempestade", \

"Nro. Vol. de Richardson", "Denom. Nro. Vol. Richardson", \
"Denom. Nro. Vol. Richardson (0-2km)", \
"Indice de Energia-Helicidade", "Parametro de Super Celula" /)

; unidade fisica dos campos
unidadeCampos = (/ "[m]","[m]","[m]","[K]","[C]","[J/kg]", \
"[J/kg]","[J/kg]","[m]",
"[m/s]","[m/s]","[m^2/s^2]","", \
"[m^2/s^2]","[m^2/s^2]","","" /)

;============= PARTE GRAFICA ========================
; criando ambiente grafico

wks = gsn_open_wks("ps","leituraBinarioETA")
; escolha do mapa de cores

gsn_define_colormap(wks,"BlWhRe")

; opcoes do grafico
Res = True
Res@gsnMaximize = True
Res@gsnAddCyclic = False

Res@gsnFrame = False
Res@gsnDraw = False

; mapa

Res@pmTickMarkDisplayMode = "Always"
Res@mpDataSetName = "Earth..4"
Res@mpDataBaseVersion = "MediumRes"

Res@mpOutlineOn = True
Res@mpOutlineSpecifiers = (/"Brazil:states"/)
Res@mpMaxLatF = max(lat)
Res@mpMinLatF = min(lat)

Res@mpMaxLonF = max(lon)
Res@mpMinLonF = min(lon)
Res@mpFillOn = False

; contorno
Res@cnFillOn = True
Res@lbLabelBarOn = True

Res@lbLabelAutoStride = True

Res@gsnSpreadColors = True
Res@cnLineThicknessF = 3.

; desenhando graficos
plot = new( 17, graphic )
do i=0,16
Res@gsnLeftString = nomeCampos(i) ; nome do campo
Res@gsnRightString = unidadeCampos(i) ; unidade fisica do campo
plot(i) = gsn_csm_contour_map( wks, saidasETA(campos|i,lat|:,lon|:), Res )

end do


; colocando tudo em duas paginas
gsn_panel( wks, plot(0:8), (/3,3/), False )

gsn_panel( wks, plot(9:16), (/3,3/), False )

; convertendo arquivo PS para arquivos PNG
system("convert -trim -density 300 -rotate -90 -geometry 1000x1000"+\
" leituraBinarioETA.ps leituraBinarioETA.png")


end


Abaixo estão os dados do arquivo binário plotados pelo script acima:

Primeiro conjunto de figuras:


Segundo conjunto de figuras:



Pronto, arquivo binário lido e informações plotadas. Gastando um pouco mais de tempo, pode-se melhorar o script para pegar as informações do CTL (dimensões das variáveis, número de variáveis, nomes e unidades físicas das variáveis etc) automaticamente.

Fica aqui a sugestão!

Até mais.

quarta-feira, 30 de setembro de 2009

Comparação NCL x GrADS: arquivos ASCII

Olá,

Abaixo mostro outra comparação simples, mas extremamente interessante.

Neste "post" coloco o que talvez seja a maior vantagem do NCL sobre o GrADS, pelo menos ao meu ver: a leitura e plotagem de dados contidos em arquivos ASCII.

Quem usa o GrADS sabe do que estou falando. Para criar um gráfico de dados contidos em um arquivo ASCII é necessário um programa (Fortran ou C) para gerar um arquivo binário dos dados do arquivo ASCII e um arquivo CTL, que possibilita a sua utilização no GrADS.

No NCL, nada disso é necessário, basta abrir o arquivo ASCII diretamente, carregando-o para uma variável (arranjo de qualquer dimensão) e plotá-lo.

Vamos aos exemplos práticos. Tenho um campo de pressão ao nível médio do mar fictício (o programa que gera este campo é apresentado no final deste "post"), armazenado em um arquivo ASCII, chamado campoFicticio.ascii. Primeiro, a tarefa de plotá-lo com o GrADS:

Passo 1: programa Fortran (90 ou 77) para criação do binário.

program campoGrads
implicit none

real, dimension(10,10) :: z

integer :: i,j

open(30, file='campoFicticio.ascii', status='old')
read(30,*) z

close(30)

open(30, file='campoFicticio.bin', &
form='unformatted', access='direct', recl=10*10*8)

write(30,rec=1) z
close(30)

end program campoGrads


Passo 2: o arquivo CTL.

DSET ^campoFicticio.bin
TITLE Campo Ficticio
UNDEF 0.10000E+16
XDEF 10 linear 180 1
YDEF 10 linear 10 1
ZDEF 1 levels 1000

TDEF 1 linear 18z30sep09 1dy

VARS 1

z 0 99 campo ficticio de pressao ao nivel do mar
ENDVARS


Passo 3: script para plotagem do campo.

'open campoFicticio.ctl'
'set grads off'
'display z'
'draw title Campo ficticio - PNMM [hPa] - GrADS'
'printim campoFicticioGRADS.png white x1000 y1000'
'quit'


O Resultado
:



Finalmente, a mesma tarefa, mas executada com o NCL.

Passo 1: leitura do arquivo ASCII e plotagem

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

begin

; lendo arquivo
arquivo = asciiread( "campoFicticio.ascii", -1, "float" )
campo = onedtond( arquivo, (/10,10/) )

; ambiente gráfico
wks = gsn_open_wks( "ps", "campoFicticio" )


; recurso gráfico, apenas para aumentar tamanho da figura
; e colocar o título
res = True

res@gsnMaximize = True

res@tiMainString = "Campo ficticio - PNMM [hPa] - NCL"

; plotando gráfico
plot = gsn_csm_contour( wks, campo, res )


; convertendo para PNG
system("convert -density 300 -geometry 1000x1000"+\
" -trim campoFicticio.ps campoFicticioNCL.png")

; apagando arquivo PS
system("rm campoFicticio.ps")

end


O Resultado:



É isso mesmo, somente um passo para plotar o mesmo gráfico no NCL. E isso serve para qualquer tipo de gráfico. Tens uma série temporal e quer plotar um gráfico XY? Mesmo processo. Basta ler a série temporal e plotá-la.

Abraços,

PS. 1: O campo usado para criar os gráficos acima foi gerado com um programa em Fortran 90, que se encontra abaixo:

program campoFicticio

! sem declaracoes implicitas
implicit none


! variaveis e parametros
integer, parameter :: nx=10, ny=10
real, parameter :: pi = acos(-1.)
real, dimension(ny,nx) :: z
integer :: i,j
real :: dx,dy

! variacao dos eixos x e y
dx = 2*pi/nx
dy = 2*pi/ny

! gerando campo ficticio
do j=0,ny-1
do i=0,nx-1
z(j+1,i+1) = sin(i*dx)*cos(j*dy)*20.+1010.
end do
end do

! escreve campo no arquivo ASCII
open(30, file='campoFicticio.ascii')
write(30,*) z
close(30)

end program campoFicticio

PS. 2: Abaixo as linhas de comando para compilar os programas f90 e rodar os scripts GrADS e NCL

FORTRAN:
$ gfortran -o programaFORTRAN.exe programaFORTRAN.f90

GRADS:
$ grads -bpc scriptGRADS.gs

NCL:
$ ncl scriptNCL.ncl

nas quais programaFORTRAN, scriptGRADS e scriptNCL são os nomes dos arquivos que contêm o programa em Fortran 90, o script do GrADS e o script do NCL, respectivamente, e gfortran é o compilador Fortran 90 que acompanha a maioria das distribuições linux. O símbolo $ indica o prompt da linha de comando do terminal do linux, não sendo necessária a sua digitação.

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.

terça-feira, 22 de setembro de 2009

A técnica de bootstrap

Na meteorologia (e em várias outras áreas) é comum nos depararmos com um conjunto de dados com comportamento distinto das distribuições teóricas, como a Distribuição Normal, por exemplo.

Um problema surge nestes casos: como aferir se a estimação de um parâmetro destes dados é representativa da população, ou seja, que teste de hipótese aplicar nestas situações. A transformação dos dados (mais conhecida como reexpressão dos dados) pode ser uma saída, fazendo com que os dados se aproximem de uma distribuição teórica. Por exemplo, dados de chuva são, geralmente, positivos e altamente concentrados nos valores mais baixos, facilmente descritos por distribuições exponenciais. Entretanto, alguns testes levam em conta, implicitamente, uma "normalidade" nos dados. Pode-se reexpressar estes dados, levando-os a uma distribuição normal, por meio da transformação logarítimica [Wilks 2006, pág. 47].

O analista pode optar, o que normalmente ocorre, por trabalhar com seus dados originais, ou seja, sem transformação alguma. É neste momento que muito cuidado deve ser tomado ao se aferir a representatividade de uma estimação. A aplicação de testes a dados que não seguem a distribuição assumida pelos primeiros pode levar a erros na rejeição ou aceitação da hipótese nula (Erros Tipo I e II, ver Wilks [2006], pág. 133).

Um alternativa aos testes paramétricos (aqueles que levam em consideração algum tipo de distribuição) é o teste não paramétrico. Há vários testes não paramétricos conhecidos, como o teste de Mann-Kendall, usado para identificação de tendências em uma série de dados e o teste de Mann-Whitney, usada para testar se duas amostras foram retiradas de populações com médias iguais, entre vários outros.

Com a maior facilidade tecnológica dos dias de hoje, uma outra abordagem ao teste não paramétrico tem ganhado cada vez mais espaço: os testes de reamostragem, que inclui o Bootstrap. A técnica do Bootstrap realiza uma reamostragem dos dados, com substituição, com o objetivo principal de obter uma distribuição do parâmetro em teste. Para maiores informações a respeito desta técnica, veja o Capítulo 5, Seção 3.4 de Wilks [2006].

O objetivo deste "post" é mostrar como usar o Bootstrap de uma maneira bem simples, usando o R (http://www.r-project.org/), um software gratuito destinado a aplicações estatísticas. Veja o script abaixo:

-------------------- INÍCIO DO SCRIPT ---------------------

# dados de precipitacao (em pol)
p <- c(0.44,1.18,2.69,2.08,3.66,1.72,2.82,0.72,1.46,1.30, 1.35,0.54,2.74,1.13,2.50,1.72,2.27,2.82,1.98,2.44,2.53,2, 1.12,2.13,1.36,4.9,2.94,1.75,1.69,1.88,1.31,1.76,2.17,2.38, 1.16,1.39,1.36,1.03,1.11,1.35,1.44,1.84,1.69,3,1.36,6.37, 4.55, 0.52,0.87,1.51)

teste_sign <- 0.025 # alpha/2, nivel de significancia.
iBoot <- 10000 # numero de reamostragens
p_log <- log(p) # logaritmo da chuva
sd_p_log <- sd(p_log) # desvio padrao dos dados originais
sd_precip_boot <- array(0, c(iBoot)) # desvios padroes reamostras

# loop que faz a reamostragem e calcula o desvio padrao dela
for (i in 1:iBoot) {
dummy <- sample(p_log,replace=T) # reamost. c/ subst.
sd_precip_boot[i] <- sd(dummy) # SD da reamostragem
rm(dummy) # apaga reamostragem
}

# calula quantis em funcao do nivel do teste de significancia
# definido na variavel 'teste_sign'
quantis <- quantile(sd_precip_boot,prob=seq(0.,1.,teste_sign))
nQuantis <- NROW(quantis) # nro de quantis calculados

# plota distribuicao (histograma) dos desvios padroes calculados
# das reamostragens
png("testeBootWilks.png",width=600,height=800) # cria PNG

# desenha histograma com 100 classes
hist(sd_precip_boot, main="Exemplo 5.10 de Wilks [2006]", xlab="desvio padrão", ylab="frequência", breaks=100)
abline(v=quantis[2],lty=2,lwd=2) # lim. inf. interv. confianca
abline(v=quantis[nQuantis-1],lty=2,lwd=2) # lim. sup. interv. confianca abline(v=sd_p_log,lty=1,lwd=2) # desvio padrao dados originais

dev.off() # fecha grafico

---------------- FIM DO SCRIPT ------------------------
As linhas tracejadas dão o Intervalo de Confiança obtido com o Bootstrap dos dados de chuva. Entre estas duas linhas está uma linha cheia, mostrando o valor do desvio padrão do logaritmo das precipitações da amostra original. Resultado: temos 95% de confiança de que o desvio padrão do logaritmo das precipitações está entre 0,41 e 0,65.

Numa pesquisa rápida nas funções disponíveis no R é possível encontrar funções que fazem diretamente o Bootstrap. Todavia, preferi mostrá-lo desta forma por ser mais didático, pois permite um entendimento melhor do processo.

Este método também pode ser aplicado a duas amostras. Mostrarei esta aplicação num "post" posterior.

Abraços,

REFERÊNCIAS:

Wilks (2006), Statistical Methods in the Atmospheric Sciences, Academic Press, 648 pp.

LEITURAS RECOMENDADAS:

LÚCIO, P. S., I. V. LEANDRO e T. P. Paula (2006), Bootstrap aplicado à avaliação de incertezas estatísticas no prognósticos de quantis extremos de precipitação, XIV Congresso Brasileiro de Meteorologia, Florianópolis, SC. Disponível em: http://www.criatividadecoletiva.net/cbm-files/14-eee1d97e2f1534538b995b911f9082b6.pdf