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.

Nenhum comentário:

Postar um comentário