Gráficos no IDL e comparação de dados

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:

Temperatura_do_ar_2013-11-12

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:

Histo_Temperatura_do_ar

Veja também o primeiro post sobre IDL.

One comment

Leave a Reply

O seu endereço de e-mail não será publicado. Campos obrigatórios são marcados com *

Esse site utiliza o Akismet para reduzir spam. Aprenda como seus dados de comentários são processados.