O IDL (sigla de Interactive Data Language) é uma linguagem de programação cientifica utilizada nas mais diversas áreas de pesquisa, conforme falado no post anterior sobre IDL. A saída gráfica do programa utiliza janelas (“windows”) para exibir seu conteúdo, cada uma identificada por um número. Seu tamanho e posição podem ser especificados, conforme o exemplo a seguir:
IDL> window, 0, xsize = 300, ysize = 400, xpos = 300, ypos = 400
Para fazer o gráfico mesmo sem criar a “window”, deve ser utilizado o “Z graphics buffer”, especificando a resolução utilizando o comando “device”:
thisDevice = !D.Name Set_Plot, 'Z' Erase Device, Set_Resolution=[300,400], Z_Buffer=0
PLOT é a rotina que irá plotar dois vetores de mesmo tamanho entre si (ou somente um vetor, com o eixo x sendo uma sequência numerada). Os parâmetros incluem: title, xtitle e ytitle (títulos), xrange e yrange (determina valores máximos e força escala), psym (símbolo), color e outros. Veja mais propriedades de gráficos no site da Exelis.
Também é possível plotar gráficos de superfície (procedimento SURFACE, parâmetros az e ax para rotacionar), com linhas ou com preenchimento (SHADE_SURF), e gráficos de contorno (procedimento CONTOUR, opções xstyle, ystyle e nlevels, que fixa o número de níveis, e palavra chave FOLLOW para imprimir os valores do contorno).
Comparação de dados
Suponha a seguinte tarefa: comparar dados de uma mesma variável, obtida nas mesmas condições, mas por diferentes instrumentos. Uma sugestão é começar fazendo séries temporais dos valores obtidos pelo instrumento de referência e do instrumento a ser comparado em um mesmo gráfico. Desse modo, pode-se visualizar a disposição temporal dos valores obtidos. Posteriormente, fazer outro gráfico com os dados do sensor teste versus sensor de referência. Outra opção é calcular as diferenças para os dados medidos simultaneamente e construir um histograma dessas diferenças. Desse modo, é possível visualizar a distribuição de frequência dessas diferenças e determinar o bias (ou viés), que é a tendência que os valores obtidos pelo instrumento assumem quando comparados ao instrumento de referência (como superestimar ou subestimar, por exemplo).
Segue um exemplo para plotar diferentes series temporais no mesmo gráfico para comparação. Essa rotina funciona para um instrumento de referência e dois instrumentos a serem comparados, onde cada equipamento é uma sequência no gráfico (Vermelho: referência, Verde: sensor1 e Preto: sensor2). Os dados de entrada devem estar separados em arquivos diários, com o nome no formato “nome_do_sensor_anomesdia.dat”, com médias a cada cinco minutos. Mesmo que dados não correspondam aos mesmos instantes temporais, ou seja, tenham falhas, o script somente plota os pontos que estejam nos arquivos.
São produzidos um gráfico para cada dia e para cada uma das seguintes variáveis: Temperatura do ar, Umidade Relativa, Pressão, Precipitação, Precipitação acumulada, Velocidade média do vento, Velocidade máxima do vento, Direção média do vento, Direção máxima do vento. Nem todas as variáveis são obtidas por todos os equipamentos.
@converte
dir_sensores = ['arq_diarios/files_sensor1/','arq_diarios/files_sensor2/']
dir_referencia = 'arq_diarios/files_referencia/'
for ano=2013,2014 do begin
for mes=1,12 do begin
for dia=1,31 do begin
converte,mes,mm$
converte,dia,dd$
yymmdd$=strmid(ano,4,4) + '-' + mm$ + '-' + dd$
flag = intarr(3)
result=strarr(2)
result_referencia = findfile(dir_referencia + '*' + yymmdd$ + '.dat',count=preferencia)
for i=0,1 do begin
result(i) = findfile(dir_sensores(i) + '*' + yymmdd$ + '.dat',count=p1)
flag(i) = p1
endfor
;Comparar sensores com referencia, caso tenha dados dele
if total(flag) gt 0 and preferencia eq 1 then begin
print,result
nlines = intarr(3)
nome_var = strarr(10)
nome_var = ['Tempo','Temperatura_do_ar','UR','Pressao','Precipitacao','Precipitacao_acumulada','Velocidade_media','Velocidade_maxima','Direcao_media','Direcao_maxima']
;Montar vetor com dados referencia
CLOSE,1
OPENR,1,result_referencia
nlines(0) = FILE_LINES(result_referencia)
julian = 0.0d
referencia = FLTARR(10,nlines(0))
serie_tempox = DBLARR(nlines(0))
i=0
while eof(1) eq 0 do begin
;2013-11-01 20:30,20.59,82.42,931.4610,0,0,1.10,8.94,104.17,104.12
;Tempo | Temperatura | UR | Pressão | Precipitação | Precipitação acumulada | Vel média | vel máx | Direção média | Direção máxima
READF,1, format='(i4,4(1x,i2),9(1x,f0))',a,m,d,hh,min,tar,ur,p,ppt,pptacum,velmed,velmax,dirmed,dirmax
;converte data pra JULDAY+fracao_do_dia
serie_tempox(i) = JULDAY(m,d,a,hh,min,0)
referencia(0,i)=serie_tempox(i)
referencia(1,i)=tar
referencia(2,i)=ur
referencia(3,i)=p
referencia(4,i)=ppt
referencia(5,i)=pptacum
referencia(6,i)=velmed
referencia(7,i)=velmax
referencia(8,i)=dirmed
referencia(9,i)=dirmax
i=i+1
endwhile
CLOSE,1
;print,referencia
;print,'----------------------------------------------------------'
;Montar vetor com dados sensor1
CLOSE,2
OPENR,2,result(0)
nlines(1) = FILE_LINES(result(0))
sensor1 = FLTARR(10,nlines(1))
serie_tempoh = DBLARR(nlines(1))
i=0
while eof(2) eq 0 do begin
;2013-11-09 00:00,18.16904,76.0881,0,NA
;Tempo | Temperatura | UR | Precipitação | Precipitação acumulada
READF,2, format='(i4,4(1x,i2),4(1x,f0))',a,m,d,hh,min,tar,ur,ppt,pptacum
;print,a,m,d,hh,min,tar,ur,ppt,pptacum
;converte data pra JULDAY+fracao_do_dia
serie_tempoh(i) = JULDAY(m,d,a,hh,min,0)
sensor1(0,i)=serie_tempoh(i)
sensor1(1,i)=tar
sensor1(2,i)=ur
sensor1(3,i)=p
sensor1(4,i)=ppt
sensor1(5,i)=pptacum
sensor1(6,i)=0
sensor1(7,i)=0
sensor1(8,i)=0
sensor1(9,i)=0
i=i+1
endwhile
CLOSE,2
;print,sensor1
;print,'----------------------------------------------------------'
;Montar vetor com dados sensor2
CLOSE,3
OPENR,3,result(1)
nlines(2) = FILE_LINES(result(1))
sensor2 = FLTARR(10,nlines(2))
serie_tempo = DBLARR(nlines(2))
i=0
while eof(3) eq 0 do begin
;2013-11-01 14:25,21.8,64.5,934,NA,2.1,5,163,243
;Tempo | Temperatura | UR | Pressão | Precipitação | Vel média | Vel máx | Direção média | Direção máx
READF,3, format='(i4,4(1x,i2),8(1x,f0))',a,m,d,hh,min,tar,ur,p,ppt,velmed,velmax,dirmed,dirmax
;print,a,m,d,hh,min,tar,ur,p,ppt,velmed,velmax,dirmed,dirmax
;converte data pra JULDAY+fracao_do_dia
serie_tempo(i) = JULDAY(m,d,a,hh,min,0)
sensor2(0,i)=serie_tempo(i)
sensor2(1,i)=tar
sensor2(2,i)=ur
sensor2(3,i)=p
sensor2(4,i)=ppt
sensor2(5,i)=0
sensor2(6,i)=velmed
sensor2(7,i)=velmax
sensor2(8,i)=dirmed
sensor2(9,i)=dirmax
i=i+1
endwhile
CLOSE,3
;print,serie_tempo, format='(d14.6)'
;print,'----------------------------------------------------------'
;Plotar grafico diario: 3 series x tempo para cada variavel
;Limites eixo y
ymin = [0,10,0,920,0,0,0,0,0,0]
ymax = [0,40,100,950,5,1,30,50,400,400]
units = ['tempo',' (oC)',' (%)',' (hPa)',' (mm)',' (mm)',' (m/s)',' (m/s)',' (o)',' (o)']
for i=1,9 do begin
serie_referencia = referencia(i,*)
serie_sensor1 = sensor1(i,*)
serie_sensor2 = sensor2(i,*)
device,decomposed=0
window,0,XSIZE=800,YSIZE=300
tvlct,0,255,0, 1 ;cor 1 = verde
tvlct,255,0,0, 2 ;cor 2 = vermelho
tvlct,0,0,255, 3 ;cor 3 = azul
date_label = LABEL_DATE( DATE_FORMAT=['%M-%D!C%H:%I'] )
plot,serie_tempo,serie_sensor2,$
XTICKFORMAT='LABEL_DATE',XTICKUNITS='Time',title=nome_var(i)+' '+yymmdd$,$
XCHARSIZE=1.3,YCHARSIZE=1.3,YTITLE=nome_var(i)+units(i),$
BACKGROUND=255,COLOR=0,LINESTYLE = 1, THICK = 2, psym = 2, YRANGE = [ymin(i), ymax(i)]
oplot,serie_tempoh,serie_sensor1,LINESTYLE = 0, THICK = 2, color = 1, psym = 1 ;cruz
oplot,serie_tempox,serie_referencia,LINESTYLE = 5, THICK = 2, color = 2, psym = 4 ;diamond
write_png,nome_var(i)+'_'+yymmdd$+'.png',tvrd(true=1)
;stop
endfor
endif
endfor
endfor
endfor
end
Observações:
– O primeiro gráfico plotado define a grade e os limites, então assume-se que o sensor cujos dados possuem a série temporal “serie_tempo” seja o mais completa o possível.
– O comando tvlct define as cores a partir do RGB, e oplot permite plotar mais de uma sequência de dados em um mesmo gráfico.
– A rotina “converte.pro”, chamada no início do script, é utilizada para incluir um zero antes do número (transforma mês 5 em 05, por exemplo), conforme segue abaixo:
pro converte,dia,dd$ if dia lt 10 then begin dd$ = '0' + strmid(fix(dia),7,1) endif else begin dd$ = strmid(fix(dia),6,2) endelse return end
A seguir, um exemplo de gráfico impresso pela rotina apresentada acima:
Muito do algoritmo utilizado no exemplo anterior é utilizado para montar os histogramas. São construídos um gráfico para cada variável, utilizando toda a série temporal.
@converte
;histograma das diferenças entre as medidas dos sensores em instantes coincidentes
nome_var = strarr(10)
nome_var = ['Tempo','Temperatura_do_ar','UR','Pressao','Precipitacao','Precipitacao_acumulada',$
'Velocidade_media','Velocidade_maxima','Direcao_media','Direcao_maxima']
units = ['tempo',' (oC)',' (%)',' (hPa)',' (mm)',' (mm)',' (m/s)',' (m/s)',' (o)',' (o)']
dir_sensores = ['arq_diarios/files_sensor1/','arq_diarios/files_sensor2/']
dir_referencia = 'arq_diarios/files_referencia/'
;Loop envolvendo variaveis
for var=1,9 do begin
nlines = 100000
dif_r1 = DBLARR(2,nlines)
dif_r2 = DBLARR(2,nlines)
dif_12 = DBLARR(2,nlines)
conth=0
contw=0
conthw=0
for ano=2013,2014 do begin
for mes=1,12 do begin
for dia=1,31 do begin
converte,mes,mm$
converte,dia,dd$
yymmdd$=strmid(ano,4,4) + '-' + mm$ + '-' + dd$
flag = intarr(3)
result=strarr(2)
result_referencia = findfile(dir_referencia + '*' + yymmdd$ + '.dat',count=preferencia)
for i=0,1 do begin
result(i) = findfile(dir_vaisala(i) + '*' + yymmdd$ + '.dat',count=p1)
flag(i) = p1
endfor
;Comparar sensores com referencia, caso tenha dados dele
if total(flag) gt 0 and preferencia eq 1 then begin
print,result
nlines = intarr(3)
;Montar vetor com dados referencia
CLOSE,1
OPENR,1,result_referencia
nlines(0) = FILE_LINES(result_referencia)
julian = 0.0d
referencia = FLTARR(10,nlines(0))
serie_tempox = DBLARR(nlines(0))
i=0
while eof(1) eq 0 do begin
;2013-11-01 20:30,20.59,82.42,931.4610,0,0,1.10,8.94,104.17,104.12
;Tempo | Temperatura | UR | Pressão | Precipitação | Precipitação acumulada | Vel média | vel máx | Direção média | Direção máxima
READF,1, format='(i4,4(1x,i2),9(1x,f0))',a,m,d,hh,min,tar,ur,p,ppt,pptacum,velmed,velmax,dirmed,dirmax
;converte data pra JULDAY+fracao_do_dia
serie_tempox(i) = JULDAY(m,d,a,hh,min,0)
referencia(0,i)=serie_tempox(i)
referencia(1,i)=tar
referencia(2,i)=ur
referencia(3,i)=p
referencia(4,i)=ppt
referencia(5,i)=pptacum
referencia(6,i)=velmed
referencia(7,i)=velmax
referencia(8,i)=dirmed
referencia(9,i)=dirmax
i=i+1
endwhile
CLOSE,1
;Montar vetor com dados sensor1
CLOSE,2
OPENR,2,result(0)
nlines(1) = FILE_LINES(result(0))
sensor1 = FLTARR(10,nlines(1))
serie_tempoh = DBLARR(nlines(1))
i=0
while eof(2) eq 0 do begin
;2013-11-09 00:00,18.16904,76.0881,0,NA
;Tempo | Temperatura | UR | Precipitação | Precipitação acumulada
READF,2, format='(i4,4(1x,i2),4(1x,f0))',a,m,d,hh,min,tar,ur,ppt,pptacum
;print,a,m,d,hh,min,tar,ur,ppt,pptacum
;converte data pra JULDAY+fracao_do_dia
serie_tempoh(i) = JULDAY(m,d,a,hh,min,0)
sensor1(0,i)=serie_tempoh(i)
sensor1(1,i)=tar
sensor1(2,i)=ur
sensor1(3,i)=p
sensor1(4,i)=ppt
sensor1(5,i)=pptacum
sensor1(6,i)=0
sensor1(7,i)=0
sensor1(8,i)=0
sensor1(9,i)=0
i=i+1
endwhile
CLOSE,2
;Montar vetor com dados sensor2
CLOSE,3
OPENR,3,result(1)
nlines(2) = FILE_LINES(result(1))
sensor2 = FLTARR(10,nlines(2))
serie_tempow = DBLARR(nlines(2))
i=0
while eof(3) eq 0 do begin
;2013-11-01 14:25,21.8,64.5,934,NA,2.1,5,163,243
;Tempo | Temperatura | UR | Pressão | Precipitação | Vel média | Vel máx | Direção média | Direção máx
READF,3, format='(i4,4(1x,i2),8(1x,f0))',a,m,d,hh,min,tar,ur,p,ppt,velmed,velmax,dirmed,dirmax
;print,a,m,d,hh,min,tar,ur,p,ppt,velmed,velmax,dirmed,dirmax
;converte data pra JULDAY+fracao_do_dia
serie_tempow(i) = JULDAY(m,d,a,hh,min,0)
sensor2(0,i)=serie_tempow(i)
sensor2(1,i)=tar
sensor2(2,i)=ur
sensor2(3,i)=p
sensor2(4,i)=ppt
sensor2(5,i)=0
sensor2(6,i)=velmed
sensor2(7,i)=velmax
sensor2(8,i)=dirmed
sensor2(9,i)=dirmax
i=i+1
endwhile
CLOSE,3
serie_referencia = referencia(var,*)
serie_sensor1 = sensor1(var,*)
serie_sensor2 = sensor2(var,*)
for k=0,nlines(1)-1 do begin
pos = where(serie_tempox eq serie_tempoh(k),npts)
if npts gt 0 then begin
j = pos
dif_r1(0,conth)=serie_tempox(j)
dif_r1(1,conth)=referencia(var,j)-sensor1(var,k)
conth = conth + 1
endif
endfor
;Busca valores do sensor2 para cada j valor da referencia
for k=0,nlines(2)-1 do begin
pos = where(serie_tempox eq serie_tempow(k),npts)
if npts gt 0 then begin
j = pos
dif_r2(0,contw)=serie_tempox(j)
dif_r2(1,contw)=referencia(var,j)-sensor2(var,k)
contw=contw+1
endif
endfor
;Busca valores do sensor1 e sensor2 para cada j valor da referencia
for k=0,nlines(0)-1 do begin
posw = where(serie_tempow eq serie_tempox(k),nptsw)
posh = where(serie_tempoh eq serie_tempox(k),nptsh)
if npts gt 0 then begin
dif_12(0,conthw)=serie_tempox(k)
dif_12(1,conthw)=sensor1(var,posh)-sensor2(var,posw)
conthw=conthw+1
endif
endfor
;Termina comparação sensores com referencia, caso tenha dados dele
endif
endfor
endfor
endfor
dif_r1 = dif_r1(*,0:conth-1)
dif_r2 = dif_r2(*,0:contw-1)
dif_12 = dif_12(*,0:contw-1)
ymin = [0,-2,-10,-10,-1,-1,-20,-20,-100,-100]
ymax = [0,+2,+10,+10,+1,+1,+20,+20,+100,+100]
bvar = [0,0.1,0.5,0.5,0.01,0.01,0.5,0.5,1,1]
units = ['tempo',' (oC)',' (%)',' (hPa)',' (mm)',' (mm)',' (m/s)',' (m/s)',' (o)',' (o)']
npts = fix((ymax(var)-ymin(var))/bvar(var)) + 1
xdif = findgen(npts)*bvar(var) + ymin(var)
histo_r1 = histogram(dif_xh,min=ymin(var),max=ymax(var),bin=bvar(var))
histo_r1 = histogram(dif_xw,min=ymin(var),max=ymax(var),bin=bvar(var))
histo_12 = histogram(dif_hw,min=ymin(var),max=ymax(var),bin=bvar(var))
device,decomposed=0
window,0,XSIZE=800,YSIZE=600
!p.multi=[0,1,3]
plot,xdif,histo_r1,xtitle=nome_var(var)+units(var),title='Referencia-Sensor1',CHARSIZE=2,BACKGROUND=255,COLOR=0
plot,xdif,histo_r2,xtitle=nome_var(var)+units(var),title='Referencia-Sensor2',CHARSIZE=2,BACKGROUND=255,COLOR=0
plot,xdif,histo_12,xtitle=nome_var(var)+units(var),title='Sensor1-Sensor2',CHARSIZE=2,BACKGROUND=255,COLOR=0
write_png,'Histo_'+nome_var(var)+'.png',tvrd(true=1)
endfor
end
A seguir, um exemplo de gráfico impresso pela rotina apresentada acima:
Veja também o primeiro post sobre IDL.
One comment