Gostaria de ajudar?
Mostrando postagens com marcador Script. Mostrar todas as postagens
Mostrando postagens com marcador Script. Mostrar todas as postagens

segunda-feira, 24 de agosto de 2015

GrADS Preenchimento

Olá a todos,

Para usuário do GrADS-2.1 é possível plotar um gráfico em que um determinado intervalo de dados pode ser destacado (hachurado) e isto se tornou possível graças a biblioteca gráfica chamada de Cairo.

A utilização desta biblioteca está bem documentada em alguns sites na internet, mas a pedidos segue um dica simples de como utilizar o comando 'set tile' para destacar um determinado intervalo dos dados. Isto pode ser útil para quem deseja destacar os dados estatisticamente significativo.

A ideia do comando "tile" é desenhar pequenas imagens para cobrir uma determinada área. O GrADS na versão 2.1 suporta a utilização do comando "tile" automaticamente, versões anteriores não possuem suporte.

Segue a descrição do comando:

'set tile nu tipo largura altura llargura fcor bcor'

onde:

nu é o número do tile que pode variar de 0 até 2000;

tipo é o tipo do tile que você quer utilizar; pode variar de 0 a 8 e ao escolher o 0 significa que você irá utilizar um arquivo com os tile (externo);

1 - sólido ;
2 - pontos ;
3 - linhas diagonais esquerda \\\\\\ ;
4 - linhas diagonais direita ////// ;
5 - linhas diagonais esquerda sobrepostas diagonais direita;
6 - linhas verticais ;
7 - linhas horizontais ;
8 - linhas horizontais sobrepostas as linhas verticais.

largura é a largura do "tile" gerada em pixels. O padrão é 9.
altura é a altura do "tile" gerada em pixels. O padrão é 9.
llargura é a espessura da linha utilizada para desenhar a linha especificada na imagem.
fcor é a cor usada para desenhar a linha
bcor é a cor do fundo utilizada para encher o fundo da imagem.

Para os "tile" gerados, uma pequena imagem é criada usando a largura e altura especificada, e, em seguida, ou é preenchido com uma cor sólida, ou é preenchido com a cor de fundo e, em seguida, uma ou duas linhas são desenhadas na imagem usando a cor do primeiro plano e a espessura da linha especificada. Note que as linhas são desenhadas de canto a canto (para linhas diagonais) ou no centro da imagem (diagonal / ou diagonal \ ). Assim, o tamanho da imagem, controlado pelos comando de largura e altura, vai controlar o espaçamento das linhas para o padrão final, e a relação de aspecto da imagem vai controlar o ângulo das linhas. Uma imagem pequena resultará em linhas espaçadas de forma mais estreita no padrão final, e uma imagem maior irá resultar num padrão mais espaçados.

Note que a espessura padrão é 3, a cor de primeiro plano padrão é o número de cores 1, e a cor de fundo padrão é um preto totalmente transparente (vermelho=0, verde=0, azul=0, alfa=0).

O "tile" será descrito como uma combinação do comando 'set tile' e 'set rgb' para definir as cores do primeiro plano e do fundo. Eles podem ser definidos em qualquer ordem, mas deve ser definido antes de qualquer desenho (display) feito na página atual.

Para utilizar um padrão de cores (rgb) ao "tile" que foi criado o comando é o seguinte:

'set rgb cc tile nu'

onde:

cc é o número do RGB criado;
nu é o número do tile.

Ao utilizar um padrão de cores a um "tile" criado não poderá ser utilizado a outro "tile".

Exemplo de utilização:

'reinit'
'sdfopen uwnd.mon.mean.nc'
'set display color white';'c'
'set grads off';'set grid off'

i = 60
n = 60
cmd = 'set rbcols'
while (i<256) 
  'set rgb 'n' 'i' 'i' 'i
  cmd = cmd%' '%n
  i = i + 15
  n = n + 1
endwhile
cmd

* definindo duas cores
'set rgb 20 255 0 0'
'set rgb 21 0 0 255'

* definindo dois "tile"
'set tile 0 2 6 6 3 20'
'set tile 1 2 6 6 3 21'
'set rgb 22 tile 0'
'set rgb 23 tile 1'

'set gxout shaded'
'd uwnd'

* Desenhando pontilhados azuis e vermelhos para valores de vento superiores a 5 e inferiores a 5.
'set clevs 5'
'set ccols -1 22'
'd uwnd'
'set clevs -5'
'set ccols 23 -1'
'd uwnd'



Exemplo e:

'reinit'
'sdfopen uwnd.mon.mean.nc'
'set display color white';'c'
'set grads off';'set grid off'
'set tile 20 8 5 5 2 0'
'set rgb 30 tile 20'
'set gxout shaded'
'd uwnd'
'cbar'
'set clevs 5'
'set ccols -1 30'
'd uwnd'
'set clevs -5'
'set ccols 30 -1'
'd uwnd'


Abaixo segue um exemplo que mostra uma série de padrões diferentes. As especificações do comando 'tile' que foi usado para criar o padrão é mostrado em cada quadrado (sem os valores de cor). As definições de cor não são mostradas. Os comandos usados para desenhar o último quadrado são mostrados abaixo. Os comandos utilizados para os outros quadrados são semelhantes.

'set rgb 36 255 255 0'
'set rgb 40 255 0 255'
'set tile 21 2 9 9 3 40 36 '
'set rgb 31 tile 21'
'set line 31'
'draw recf 7 1 9 3 '
'draw string 8 0.8 set tile 21 2 9 9 3'



Boa Sorte e Bom Trabalho a Todos!



segunda-feira, 22 de junho de 2015

Plot TSM, Gelo Marinho e Cobertura de Neve - GrADS


Você dificilmente vai fazer uma determinada figura se você não precisa dela. Recentemente eu precisei plotar a cobertura de Gelo Marinho a partir de dados de TSM. Vasculhando a internet encontrei um exemplo de um plot bem legal:

Os dados utilizados neste exemplo são do GFS e para ler os dados mais recentes eu irei utilizar o exemplo mostrado aqui no blog no post anterior.

Para que este script funcione perfeitamente é necessário ter as seguinte funções do Kodama em sua maquina: color.gs, xcbar.gs e colormap.gs.

Os primeiros comandos irão pegar a informação de data no sistema e abrir os dados de TSM e Gelo do GFS com resolução 0.25.

'reinit'
'! date +%Y%m%d > date.txt'
date=read('date.txt')
vardate=sublin(date,2)

arqgfs='http://nomads.ncep.noaa.gov:9090/dods/gfs_0p25/gfs'vardate'/gfs_0p25_00z'
arqgelo='http://nomads.ncep.noaa.gov:9090/dods/ice/ice'vardate'/ice.00z'

'sdfopen 'arqgfs
'sdfopen 'arqgelo

'!rm date.txt'

Com os dados lidos iremos "setar" algumas funções para plotar os dados utilizando a projeção robinson.

'set gxout shaded'
'set mpdset hires'
'set lon -180 180'
'set lat -90 90'
'set grads off'
'set map 0 1 6'
'set mproj robinson'
'colormaps -l 272 307 0.5 -map jet'
'd tmpsfc'
'xcbar -fs 4'

Após estes comandos o plot exibido é +/- assim:



Agora iremos criar uma mascara para os continentes:

'set rgb 73 80 80 80'
'basemap L 73 0 M'

O que geramos a partir destes comandos é algo assim:



Em seguida iremos plotar o Gelo Marinho, para isto utilizamos os comandos:

'set gxout shaded'
'set dfile 2'
'set z 1'
'set t 1'
'set map 0 1 6'
'color 0 0.6 0.1 -kind dimgray->seashell->white'
'd maskout(icecmsl,icecmsl-0.05)'

O que iremos obter é o seguinte plot:



Por último iremos plotar a Cobertura de Neve e os últimos comandos do script são os seguintes:

'set dfile 1'
'color 5 100 5 -kind dimgray->seashell->white'
'set z 1'
'set t 1'
'set map 0 1 6'
'd maskout(weasdsfc,weasdsfc-5)'

E a última figura exibida tem o seguinte formato:



Espero que seja útil para vocês. 

Boa Sorte e Bom trabalho a Todos!

segunda-feira, 15 de junho de 2015

Obtendo a saída do comando Shel - GrADS


Olá a Todos,

Algumas pessoas utilizam o GrADS para ler os dados do GFS direto no GrADS Data Server. Em alguns casos pode ser necessário executar um script operacionalmente para efetuar esta tarefa. Um problema que pode existir para ler este tipo de dado (GrADS Data Server) operacionalmente é que por ser um dado de modelo os nomes dos arquivos mudam constantemente e torna-se necessário criar uma variável data para atualizar o nome do dado e efetuar a leitura corretamente.

Por exemplo:

Se você deseja ler o dado do GFS 0.25 deg starting from 00Z07jun2015 a url é a seguinte:

http://nomads.ncep.noaa.gov:9090/dods/gfs_0p25/gfs20150607/gfs_0p25_00z

Observe que um dos diretórios onde está armazenado o dado contém a data (em negrito) em que o dado foi gerado ou o modelo foi rodado. Para gerar esta data automaticamente você pode criar um script SHELL e utilizar o comando date e resolver este problema, mas você pode utilizar o GrADS e resolver este problema diretamente no script .gs.

Para isso você pode utilizar o comando date no próprio script .gs . O script fica assim:

'reinit'
'! date +%Y%m%d > date.txt'
date=read('date.txt')
vardate=sublin(date,2)

*Neste caso vardate=20150607
*Observe que o comando date produz uma data diferente pois busca a data do sistema

file='http://nomads.ncep.noaa.gov:9090/dods/gfs_0p25/gfs'vardate'/gfs_0p25_18z'

'sdfopen 'file

'!rm date.txt'

Eu não conheço uma outra forma de transformar o comando SHELL em variável, a única forma que conheço é escrevendo a saída do comando SHELL em um arquivo e depois ler este arquivo em seguida criando uma variável. Após a leitura do conteúdo do arquivo date.txt eu sugiro remover ('!rm date.txt') este arquivo para não ocupar espaço.

PS: Os dados do GFS são gerados de 6 em 6 horas para criar uma variável para ler os 4 arquivos você pode gerar as variáveis e utilizar um loop (WHILE) para ler cada arquivo.

Boa Sorte e Bom Trabalho a Todos!!

segunda-feira, 8 de junho de 2015

Duas variáveis em um mesmo campo - GrADS


Olá a todos,

É importante percebermos que ao criar uma figura que será inserida em um artigo para ser publicado precisamos dedicar um pouco de esforço para fazer uma figura com a maior qualidade possível. Muitas vezes seu artigo será mais lembrado pelas figuras do que pelo conteúdo (o conteúdo é mais importante, não é isto que esta em discussão).

PS. Dica para quem tem o GrADS versão 2.

Esta dica é bastante simples, mas pode ser bastante útil quando se deseja criar figuras com alta qualidade. Existem vários comandos no GrADS que te possibilitam gerar uma figura com alta qualidade. No entanto, quando plotamos duas variáveis (sobrepostas) no mesmo campo esquecemos que no GrADS a segunda variável apesar de sobreposta é uma segunda figura (com novos eixos, grades e etc) e por isso é importante inibir a plotagem dos labels dos eixos, as linhas que representam as grades e etc desta segunda variável. Ao inibir o plot do label dos eixos por exemplo você mantêm a qualidade da font das Latitudes e Longitudes por exemplo. Para que isso não ocorra use os comandos abaixo sempre antes do display da próxima variável exibida:

 'd press'
 'set xlab off'
 'set ylab off'
 'set grid off'
 'set mpdraw off'; *(Este comando inibe o plot do contorno dos continentes)
 'd temp'

 *Com o comando 'd temp' apenas os dados de temperatura serão plotados em seu campo. Espero que ajude. 

 Boa Sorte e Bom Trabalho a Todos!!

quinta-feira, 30 de abril de 2015

Imagem de Satelite no Grads

Pergunta do Arthur Lucas,

Olá Cristiano,

tenho visto em alguns artigos com campos meteorológicos sobrepostos a imagens de Satélites, é possível fazer isso no Grads ?

Obrigado,

É sim possível desde que o dado da imagem de satélite seja disponibilizado no formato que seja possível o Grads ler, por exemplo: Binário, Netcdf e etc.

Alguns Institutos disponibilizam algumas imagens em formato binário. Durante meu doutorado eu utilizei dados de temperatura de brilho e consegui plotar algumas figuras de TB com vento e etc. Segue um exemplo: Eu plotei a temperatura de brilho com o shaded do Grads mas e com um RGB em tons de cinza que fica "quase igual" a uma imagem de satélite.





Abraço e Bom Trabalho a Todos!!


sexta-feira, 3 de abril de 2015

Gerando mascara com o cdo

Ontem eu mostrei aqui como criar uma mascara (maskout) no Grads para plotar vetor vento significativo. 

O Guilherme Martins sempre que possível contribui com alguma sugestão ou script. Desta vez ele compartilha conosco uma forma de resolver este mesmo problema sendo que com o cdo (Climate Data Operator).

PS: Não deixem de olhar a Apostila do próprio Guilherme sobre cdo.

Sem mais delongas, segue:

Dado para fazer este exemplo:

u -> uwnd.nc
v -> vwnd.nc

1) Juntando os arquivos de u e v para calcular a velocidade do vento

cdo merge uwnd.nc vwnd.nc uv.nc

2) Calculando a velocidade do vento 

cdo expr,'vel=sqrt(uwnd*uwnd+vwnd*vwnd);' uv.nc vel.nc

3) Criando a máscara. Defino o intervalo da velocidade entre 0 e 4 m/s.

cdo setrtomiss,0,4 vel.nc mascara.nc

4) Aplicando a máscara na componente zonal (u) e meridional (v)

cdo ifthen mascara.nc uwnd.nc vento.u.nc

cdo ifthen mascara.nc vwnd.nc vento.v.nc

Ao plotar o que obtemos é o seguinte:


Boa Sorte e Bom Trabalho a Todos!

quinta-feira, 2 de abril de 2015

Grads Plot do Vento Significante


Recentemente precisei plotar somente os vetores de vento no Grads significativos, ou seja, dentro de uma área/região eu necessitava que apenas um determinado limiar de magnitude do vetor vento fosse exibido.

Uma solução para este problema é utilizar uma máscara (maskout) para "mascarar" um determinado limiar de magnitude que você não quer exibir. Portanto, o script ou parte dele fica mais ou menos assim:

'open /usr/local/opengrads/Resources/SampleDatasets/model.ctl'

'set display color white'
'c'
'set lev 850'
'set t 1'
'set lon -100 50'
'set lat 0 50'

'wind = mag(ua,va)'

** Escolha do limiar que se deseja "mascarar"

limiar_vento = 8
'sigu = maskout(ua,wind-'limiar_vento')'
'sigv = maskout(va,wind-'limiar_vento')'
'sigwind = mag(sigu,sigv)'

'set gxout shaded'
'd sigwind'
'cbar'

'set gxout vector'
'd sigu;sigv'

Este seria o resultado.

Boa Sorte e Bom Trabalho a Todos!

quarta-feira, 18 de fevereiro de 2015

Alterando os níveis verticais com CDO

Sugestão do Guilherme Martins!

O CDO sempre me surpreendendo!!!!
Alterando os níveis verticais com CDO.
Os níveis verticais dos modelos do CMIP5 estão em "Pa", e normalmente utilizamos em "hPa". Vamos usar o operador "setzaxis" para realizar essa alteração.
Carreguei neste tópico o arquivo "exemplo.nc" para vocês fazerem o teste.
Informações vista com "ncdump -h exemplo.nc"
double plev(plev) ;
plev:bounds = "plev_bnds" ;
plev:units = "Pa" ; <===== unidade do nível vertical
plev:axis = "Z" ;
plev:positive = "down" ;
plev:long_name = "pressure" ;
plev:standard_name = "air_pressure" ;
Informações vista com "ncdump arquivo.nc | more"
plev = 100000, 85000, 70000, 50000, 25000, 10000, 5000, 1000 ; => "Pa"
Segue a dica:
Crie um arquivo texto, e adicione as 6 linhas abaixo:
zaxistype = pressure
size = 8
name = lev
longname = pressure
units = hPa
levels = 1000 850 700 500 250 100 50 10
Salve com um nome qualquer, por exemplo, "nivel". Esse arquivo será lido pelo operador "setzaxis":
Utilize o comando abaixo:
cdo setzaxis,nivel exemplo.nc saida.nc
Esse comando altera a estrutura do seu arquivo netCDF com as novas informações em "hPa".
Para visualizar as novas alterações:
ncdump -h saida.nc

Arquivo: exemplo.nc

PS: caso o arquivo tenha dados UNDEF (que é caso), utilize o operador fillmiss ==> cdo -fillmiss -setzaxis,nivel exemplo.nc saida.nc

Boa Sorte e Bom Trabalho a Todos!

terça-feira, 28 de outubro de 2014

cdo merge levels

Olá a Todos,

Recentemente estive tentando usar o comando merge do cdo para criar um arquivo com vários níveis. O problema era o seguinte:

Eu possuo 10 arquivos de vento em formato binário sendo cada arquivo contendo 100 tempos e um nível.

A primeira coisa que eu fiz foi criar arquivos netcdf a partir dos arquivos descritores .ctl através do comando cdo import_binary -f nc *.ctl. No entanto desta forma o cdo não definia o nível referente ao arquivo e sempre atribuía nível 0 (zero). Com os arquivos criados com a descrição de nível 0 (zero) ao juntar (merge) os 10 arquivos o resultado era sempre um arquivo com somente um nível.

Eu não consegui resolver este problema de maneira mais simples o que fiz foi o seguinte:

Utilizei o comando sdfwrite do GrADS para gerar os arquivos netcdf que precisava.

Ex:

set x 1 144
set y 1 73
set t 1 100
set z 1
define uwnd = uwnd
set sdfwrite -4d arquivo1.nc
sdfwrite uwnd
quit

OBS1: É necessário inserir a função -4d do comando sdfwrite, desta forma você força que o arquivo de saída tenha 4 dimensões (lat, lon, tempo, nível)

OBS2: O nível no arquivo ctl precisa ser definido de maneira correta, por exemplo:

          zdef 1 linear 200 1

          Qualquer definição diferente disso deu errado.

Agora que possuía os 10 arquivos .nc utilizei o merge para gerar meu arquivo com os 10 níveis e 100 tempos:

cdo merge arquivo1.nc arquivo2.nc ... arquivo10.nc saida.nc

Boa Sorte e Bom Trabalho a Todos!!

terça-feira, 7 de outubro de 2014

Dado Texto ou ASCII no GrADS

Olá Pessoal,

Através do GrADS é possível escrever e ler arquivos textos e em algum momento isto pode ser de grande utilidade. Recentemente a Giulia de Salve me escreveu com uma dúvida exatamente sobre este contexto.

O Guilherme Martins em sua Apostila já colocou um exemplo bem explicativo e também existe um guia no site do GrADS que é mais explicativo ainda. O que eu vou colocar aqui é exatamente o que esta na página neste guia do site GrADS.

PS. Eu irei utilizar como exemplo o arquivo model.ctl para quem possui o GrADS 2 (opengrads) instalado ele fica na pasta (grads-2.0.1/Contents/Resources/SampleDatasets/).

O GrADS não possui uma função intrínseca para ler arquivos texto como o sdfopen, por exemplo, por isso é necessário alguns ajustes para isso:

Imagine que você possua um arquivo texto/ascii (my_ascii_file.txt) que tenha uma coluna com valores de 1 a 30:

1
2
3
4
...

Você pode transformar este arquivo em um binário e criar um ctl para ler no GrADS.

'open model.ctl'
file = 'my_ascii_file.txt'
datafile = 'my_binary_file.dat'
'set gxout fwrite'
'set fwrite -ap 'datafile
'!/bin/rm -f 'datafile
while (1)
         res = read(file)
         line1 = sublin(res,1)
         line2 = sublin(res,2)
         rc1 = subwrd(line1,1)
         if (rc1); break; endif
              val = subwrd(line2,1)
              'd 'val
endwhile
rc = close(file)
'disable fwrite'

O ctl para ler o dado gerado (my_binary_file.dat), imaginando que você queira plotar os dados com relação ao eixo X, ficaria assim:

dset ^my_binary_file.dat
title Sample of ASCII data converted to binary
undef -9.99e8
xdef 100 linear 1 1
ydef 1 linear 1 1
zdef 1 linear 1 1
tdef 1 linear 01jan0001 1dy
vars 1
a 0 99 ascii variable
endvars

Se você abrir o ctl acima e plotar vai gerar uma figura semelhante a esta:

Uma outra opção é criar arquivos textos a partir de qualquer dado binário através da função 'set gxout print'. A formatação desta função pode ser feita através da função 'set prnopts'. A saída é plotada no terminal mas pode ser armazenada em uma variável e posteriormente através do script ser escrita em um arquivo. O exemplo a seguir escreve 72 números por linha que representa os 72 pontos de grade em X.

'open model.ctl'
'set x 1 72'
outfile = my_ascii_file.txt'
'set gxout print'
'set prnopts %7.3f 72 1'
'd ts'
rc = write(outfile,result)

Desta forma o arquivo my_ascii_file.txt irá conter algo como isso:

Printing Grid -- 3312 Values -- Undef = -9.99e+08
258.493 258.493 258.493 258.493 258.493 258.493 258.493 258.493 258.493 258.493 258.493 <etc.>

No caso anterior, observe que a primeira linha do arquivo de saída contem algumas informações sobre o número de valores escritos e valores indefinidos que você pode não querer que seja escrito. O exemplo a seguir omiti essa linha, ou seja, pula a primeira linha de informações de diagnóstico. Além disso neste exemplo a saída será delimitada por vírgulas para importar para uma planilha (excell) por fim a variável apresentada agora é uma expressão em vez de um nome. 

'open model.ctl'
'set x 1 72'
outfile='my_ascii_file.txt'
'!/bin/rm -f 'outfile
'set gxout print'
'set prnopts %g, 72 0'
'd ts-273.15'
i=1
while (1)
    line = sublin(result,i)
    if (line = ''); break; endif
    if (i>1)
    rc = write(outfile,line,append)
    endif
    i=i+1
endwhile

Neste caso o arquivo my_ascii_file.txt conterá algo como:

-14.6567,-14.6567,-14.6567,-14.6567,-14.6567,-14.6567,-14.6567,-14.6567,-14.6567,<etc.>

Obs. A função append serve para inserir os dados no final do arquivo caso ele exista. Se você não utilizar a função append e existir um arquivo de mesmo nome ele será sobrescrito. Para mais informações leia "append".

Por fim é possível escrever informações sobre a grade (lat x lon) além dos dados de mais de uma variável, para isso é necessário salvar o resultado da exibição de cada variável/expressão (incluindo longitude e latitude), para cada ponto de grade individualmente, organizar a estrutura que será escrito no arquivo texto de saída e depois escrevê-lo. Neste exemplo será escrito a longitude, latitude, temperatura da superfície, pressão de superfície e precipitação:

'open model.ctl'
'set x 1 72'
outfile='my_ascii_file.txt'
'!/bin/rm -f 'outfile
'set gxout print'
fmt='%8.3f'
numcols=72
'set prnopts 'fmt' 'numcols' 1'
'd lon'
lon_data = result
'd lat'
lat_data = result
'd ts'
v1_data = result
'd ps'
v2_data = result
'd p'
v3_data = result
i=1
while (1)
    lons = sublin(lon_data,i)
    lats = sublin(lat_data,i)
    line1 = sublin(v1_data,i)
    line2 = sublin(v2_data,i)
    line3 = sublin(v3_data,i)
    if (lons='' | lats='' | line1='' | line2='' | line3=''); break; endif
    if (i>1)
        j=1
        while (j<=numcols)
            str = subwrd(lons,j); lon = math_format(fmt,str)
            str = subwrd(lats,j); lat = math_format(fmt,str)
            str = subwrd(line1,j); v1 = math_format(fmt,str)
            str = subwrd(line2,j); v2 = math_format(fmt,str)
            str = subwrd(line3,j); v3 = math_format(fmt,str)
            record = lon' 'lat' 'v1' 'v2' 'v3
            rc = write(outfile,record,append)
        j=j+1
        endwhile
    endif
i=i+1
endwhile

Obs. A função append serve para inserir os dados no final do arquivo caso ele exista. Se você não utilizar a função append e existir um arquivo de mesmo nome ele será sobrescrito. Para mais informações leia "append".

Neste caso a saída my_ascii_file.txt apresentará a seguinte estrutura:

0.000     -90.000     258.493     669.911     0.000
5.000     -90.000     258.493     669.911     0.000
10.000   -90.000     258.493     669.911     0.000
15.000   -90.000     258.493     669.911     0.000
20.000   -90.000     258.493     669.911     0.000
25.000   -90.000     258.493     669.911     0.000
30.000   -90.000     258.493     669.911     0.000
35.000   -90.000     258.493     669.911     0.000
40.000   -90.000     258.493     669.911     0.000

A um tempo atrás coloquei a disposição uma função que ajuda a escrever um arquivo texto sem muitas complicações. No entanto o arquivo de saída não possui formatação, veja aqui.

Boa Sorte e Bom Trabalho a Todos!

Cristiano






quinta-feira, 11 de setembro de 2014

Barra de cores "pulando" intervalos

Está é mais uma sugestão de barra de cores (shaded) para usuários Grads. O interessante desta barra é que é possível "saltar" os valores (números) referentes aos contornos da variável.

O autor deste script chama-se Stephen R. McMillan


*
*  Script to plot a colorbar
*
*  The script will assume a colorbar is wanted even if there is
*  not room -- it will plot on the side or the bottom if there is
*  room in either place, otherwise it will plot along the bottom and
*  overlay labels there if any.  This can be dealt with via
*  the 'set parea' command.  In version 2 the default parea will
*  be changed, but we want to guarantee upward compatibility in
*  sub-releases.
*
*       modifications by mike fiorino 940614
*       rvp added by diane stokes 20040729
*       cinc added by diane stokes 20050130
*
*       - the extreme colors are plotted as triangles
*       - the colors are boxed in white
*       - input arguments in during a run execution:
*
*       run cbarn cinc sf vert xmid ymid rvp
*
*       cinc - increment for labels
*       sf   - scale the whole bar 1.0 = original 0.5 half the size, etc.
*       vert - 0 FORCES a horizontal bar = 1 a vertical bar
*       xmid - the x position on the virtual page the center the bar
*       ymid - the x position on the virtual page the center the bar
*       rvp  - 1 resets virtual page (one colorbar for several plots)
*
*       if cinc is not specified, every color interval is labeled
*       if vert,xmid,ymid are not specified, they are selected
*       as in the original algorithm
*
*       if rvp is not specified, it is 0 (do not reset virtual page)

function colorbar (args)

*  Parse command line arguments
*  ----------------------------
   cinc=subwrd(args,1)
   sf=subwrd(args,2)
   vert=subwrd(args,3)
   xmid=subwrd(args,4)
   ymid=subwrd(args,5)
   rvp=subwrd(args,6)
   if(cinc=''|cinc=0);cinc=1;endif
   if(sf='');sf=1.0;endif

*  Check shading information
*  -------------------------
   'query shades'
   shdinfo = result
   if (subwrd(shdinfo,1)='None')
      say 'Cannot plot color bar: No shading information'
      return
   endif

* Check if reset of virtual page is desired
* -----------------------------------------
  if(rvp = 1)
     say 'Colorbar based on colors from last plot drawn'
     'set vpage off'
  endif

* Get plot size info
* ------------------
  'query gxinfo'
  rec2 = sublin(result,2)
  rec3 = sublin(result,3)
  rec4 = sublin(result,4)
  xsiz = subwrd(rec2,4)
  ysiz = subwrd(rec2,6)
  ylo = subwrd(rec4,4)
  xhi = subwrd(rec3,6)
  xd = xsiz - xhi

  ylolim=0.6*sf
  xdlim1=1.0*sf
  xdlim2=1.5*sf
  barsf=0.8*sf
  yoffset=0.2*sf
  stroff=0.05*sf
  strxsiz=0.12*sf
  strysiz=0.13*sf

*  Decide if horizontal or vertical color bar
*  and set up constants.
*  ------------------------------------------
   if (ylo<ylolim & xd<xdlim1)
      say "Not enough room in plot for a colorbar"
      return
   endif
   cnum = subwrd(shdinfo,5)

*  Logic for setting the bar orientation with user overides
*  --------------------------------------------------------
   if (ylo<ylolim | xd>xdlim1)
       vchk = 1
       if(vert = 0) ; vchk = 0 ; endif
   else
       vchk = 0
       if(vert = 1) ; vchk = 1 ; endif
   endif

*  Vertical bar
*  ------------
   if (vchk = 1 )
     if(xmid = '') ; xmid = xhi+xd/2 ; endif
     xwid = 0.2*sf
     ywid = 0.5*sf
     xl = xmid-xwid/2
     xr = xl + xwid
     if (ywid*cnum > ysiz*barsf)
        ywid = ysiz*barsf/cnum
     endif
     if(ymid = '') ; ymid = ysiz/2 ; endif
     yb = ymid - ywid*cnum/2
     'set string 1 l 5'
     vert = 1

*  Horizontal bar
*  --------------
   else

     ywid = 0.4
     xwid = 0.8
     if(ymid = '') ; ymid = ylo/2-ywid/2 ; endif
     yt = ymid + yoffset
     yb = ymid
     if(xmid = '') ; xmid = xsiz/2 ; endif
     if (xwid*cnum > xsiz*barsf)
         xwid = xsiz*barsf/cnum
     endif
     xl = xmid - xwid*cnum/2
     'set string 1 tc 5'
     vert = 0
   endif

*  Plot colorbar
*  -------------
   'set strsiz 'strxsiz' 'strysiz
   num = 0
   while (num<cnum)
     rec = sublin(shdinfo,num+2)
     col = subwrd(rec,1)
     hi = subwrd(rec,3)
     if (vert)
        yt = yb + ywid
     else
        xr = xl + xwid
     endif

*   Draw the left/bottom triangle
*   -----------------------------
    if (num = 0)
       if(vert = 1)
          xm = (xl+xr)*0.5
          'set line 'col
          'draw polyf 'xl' 'yt' 'xm' 'yb' 'xr' 'yt' 'xl' 'yt
          'set line 1 1 5'
          'draw line 'xl' 'yt' 'xm' 'yb
          'draw line 'xm' 'yb' 'xr' 'yt
          'draw line 'xr' 'yt' 'xl' 'yt
       else
          ym = (yb+yt)*0.5
          'set line 'col
          'draw polyf 'xl' 'ym' 'xr' 'yb' 'xr' 'yt' 'xl' 'ym
          'set line 1 1 5'
          'draw line 'xl' 'ym' 'xr' 'yb
          'draw line 'xr' 'yb' 'xr' 'yt
          'draw line 'xr' 'yt' 'xl' 'ym
       endif
   endif

*  Draw the middle boxes
*  ---------------------
   if (num!=0 & num!= cnum-1)
      'set line 'col
      'draw recf 'xl' 'yb' 'xr' 'yt
      'set line 1 1 5'
      'draw rec  'xl' 'yb' 'xr' 'yt
   endif

*  Draw the right/top triangle
*  --------------------------- 
   if (num = cnum-1)
      if (vert = 1)
         'set line 'col
         'draw polyf 'xl' 'yb' 'xm' 'yt' 'xr' 'yb' 'xl' 'yb
         'set line 1 1 5'
         'draw line 'xl' 'yb' 'xm' 'yt
         'draw line 'xm' 'yt' 'xr' 'yb
         'draw line 'xr' 'yb' 'xl' 'yb
      else
         'set line 'col
         'draw polyf 'xr' 'ym' 'xl' 'yb' 'xl' 'yt' 'xr' 'ym
         'set line 1 1 5'
         'draw line 'xr' 'ym' 'xl' 'yb
         'draw line 'xl' 'yb' 'xl' 'yt
         'draw line 'xl' 'yt' 'xr' 'ym
      endif
   endif

*  Put numbers under each segment of the color key
*  but only if this level is to be printed (based on cinc)
*  -------------------------------------------------------
   if(num < cnum-1)
*     base prt on hi if increment is based on data value (less reliable)
*     prt = math_fmod(hi,cinc);
*     ------------------------------------------------------------------
      prt = math_fmod(num,cinc);
      if(prt = 0)
         if (vert)
            xp=xr+stroff
            'draw string 'xp' 'yt' 'hi
         else
            yp=yb-stroff
            'draw string 'xr' 'yp' 'hi
         endif
      endif
    endif

*   Reset variables for next loop execution
*   ---------------------------------------
    if (vert)
      yb = yt
    else
      xl = xr
    endif
    num = num + 1

  endwhile
  return

Para utilizar é fácil:

Após abrir o arquivo desejado, definir a função gxout shaded e exibir (display) a variável, basta digitar:

ga-> run cbarnskip.gs 2

Eu chamei de cbarnskip.gs você escolhe o nome que achar melhor. O número 2 está relacionado aos valores que serão "removidos" ou alternados.

O resultado pode ser algo do tipo:



Boa Sorte e Bom Trabalho a Todos!

quinta-feira, 6 de fevereiro de 2014

GxYAT um comando para gerar Figuras com alta resolução

A pedido de Luana Pampuch segue um exemplo para salvar no GrADS figuras em .pdf ou .eps e etc.

No OpenGrADS existe um comando chamado GxYAT que você pode criar figuras em PDF, PNG, Postscript e SVG.

Para quem escreve textos em Latex (meu caso) é excelente pois as figuras ficam com excelente resolução e durante a impressão as figuras não perdem qualidade.

Portanto, para quem esta escrevendo TCC, Dissertação ou Tese fica a dica, aprendam Latex salvem suas figuras em PDF e se livrem de serem cobrados pela má qualidade das figuras durante as defesas :).

como usar no GrADS:

'enable print nome_da_figura.gmf'
.
.
.
'print'
'disable print'


'quit'

No Terminal Linux digite:

> gxyat -o nome_da_figura.pdf nome_da_figura.gmf
É isso!!

Boa Sorte e Bom Trabalho a Todos

terça-feira, 4 de fevereiro de 2014

Script GrADS Barra em formato de Diamante


Olha o GrADS surpreendendo novamente. Recentemente encontrei um script que plota uma barra de cores em formato de diamante e irei compartilhar com vocês. Acredito que fazer figuras para apresentações com este tipo de barra pode ser interessante, vamos lá:


Para usar as opções são as seguintes:

'run diamondbar.gs -r -marg 0.1 -font 0.09 -fs 2'

Este comando serve para plotar a figura acima seguindo o exemplo disponibilizado abaixo.

Opções:

-help                = Ajuda;
-loc                  = Escolhe um local para a colorbar (r=direita-vertical, l=esquerda-vertical,
                           b=base-horizontal, t=topo-horizontal),
                           você também pode chamar o comando através de (-l, -r, -b, -t);
-size/s              = Escolha o tamanho dos diamantes;
-marg              = Escolha a margem para o eixo;
-fs/int              = Escolha o intervalo para cada número plotado ao longo da colorbar;
-font/fontsize   = Escolha o tamanho da fonte dos números;
-scol                 = Escolha a cor dos números;
-line                 = Desenha o contorno de cada diamante caso deseje.

Scripts utilizados:

barradiamante.gs
Exemplo.gs

Simples assim, espero que funcione.

Boa Sorte e Bom Trabalho a Todos!!

segunda-feira, 11 de novembro de 2013

Script Download Dados ECMWF/ERA-Interin - Parte 1


Quando eu atualizei o link do "novo" servidor de download de dados da Reanalise do Era-Interim o Helber comentou de um script que possibilitava o download. Segue como fazer:

PS. Para entender os parâmetros e como modificá-los não esqueça de lê o segundo artigo (Script Download Dados ECMWF/ERA-Interim - Parte 2), link no final deste post.

O ECMWF está disponibilizando uma API (Aplication Programming Interface) nas linguagens (Python, Perl e Java) para acesso/download dos seus dados via script/terminal. É uma ferramenta extremamente útil para quem esta familiarizado com este tipo de sistema e não quer baixar os dados via site. Eu recomendo!

Eu irei tentar ensinar como usar a API desenvolvida em Python, escolha a linguagem que você melhor trabalha. Todas estas informações estão disponíveis no site da API no ECMWF.

Para instalar a biblioteca (ecmwfapi) você precisará primeiro instalar o programa pip.

O programa pip é uma excelente alternativa para instalar bibliotecas desenvolvidas em python de maneira fácil, através do terminal. Para instalar o pip no Ubuntu execute no terminal os seguintes comandos:

> sudo apt-get install python-pip python-dev build-essential 
> sudo pip install --upgrade pip 

Agora com o pip instalado execute no terminal: (tudo junto separado por espaço)

> sudo pip install https://software.ecmwf.int/wiki/download/attachments/23694554/ecmwf-api-client-python.tgz

Este comando irá instalar a biblioteca, você pode verificar ela instalada no diretório:

/usr/local/lib/python2.7/dist-packages/

Feito isso, agora você irá precisar criar a sua chave de acesso. Quando precisávamos baixar os dados da reanálise do ECMWF assim que escolhíamos os dados o site solicitava os dados pessoais (Nome a Instituição e o email), agora com essa API você gera uma chave, insere no script, e não precisará fazer esta autorização sempre que for baixar, agora ela feita automaticamente via API.

Para criar a chave você precisa se cadastrar no site:

http://apps.ecmwf.int/registration/

Um email de conformação será enviado com seu login e senha, confirme clicando no link. Ao fazer login na pagina abaixo será disponibilizado sua chave de acesso:

https://apps.ecmwf.int/auth/login/

Agora que sua chave já esta disponível é preciso criar um arquivo com o nome (.ecmwfapirc) no seu /home/usuario para o acesso ser automático, Obs: o nome do arquivo é .ecmwfapirc com (.) para ficar oculto e sem extensão.

Neste arquivo você irá inserir:

{
    "url"   : "https://api.ecmwf.int/v1",
    "key"   : "XXXXXXXXXXXXXXXXXXXXXX",
    "email" : "john.smith@example.com"
}

,onde vocé irá substituir XXXXXXXXXXXXXXXXXXXXXX pela sua chave e john.smith@example.com pelo email que você utilizou para efetuar o cadastro.

Feito isso você já pode baixar os dados. O ecmwf disponibiliza um script exemplo para teste. Portanto, crie um script em PYTHON e dentro dele insira as linhas: eu chamei o meu script de exemplo.py


#!/usr/bin/env python
from ecmwfapi import ECMWFDataServer

server = ECMWFDataServer()

server.retrieve({
    'dataset' : "tigge",
    'step'    : "24/to/120/by/24",
    'number'  : "all",
    'levtype' : "sl",
    'date'    : "20071001/to/20071003",
    'time'    : "00/12",
    'origin'  : "all",
    'type'    : "pf",
    'param'   : "tp",
    'area'    : "70/-130/30/-60",
    'grid'    : "2/2",
    'target'  : "data.grib"
    })

Salve e dê permissão de execução:

> chmod +x exemplo.py 

Ao executar este script:

>./exemplo.py

Aparecerá uma série de informações no terminal sobre o sistema, dado e o andamento do Download.

Script Download Dados ECMWF/ERA-Interim - Parte 2

Boa Sorte e Bom Trabalho a Todos.

Script Download Dados ECMWF/ERA-Interin - Parte 2

Agora que já sabemos como instalar a API do ECMWF e já fizemos o teste do download vamos explicar o que cada parâmetro do script (exemplo.py) significa e principalmente se eu precisar baixar um dado onde devo modificar, o que modificar?

No script (exemplo.py) do post anterior:

#!/usr/bin/env python
from ecmwfapi import ECMWFDataServer

server = ECMWFDataServer()

server.retrieve({
    'dataset' : "tigge",
    'step'    : "24/to/120/by/24",
    'number'  : "all",
    'levtype' : "sl",
    'date'    : "20071001/to/20071003",
    'time'    : "00/12",
    'origin'  : "all",
    'type'    : "pf",
    'param'   : "tp",
    'area'    : "70/-130/30/-60",
    'grid'    : "2/2",
    'target'  : "data.grib"
    })

dataset = servidor de dado. O ERA-interim se chama interim (independente se for superfície ou níveis de pressão);
step = time step. Neste exemplo tigge são dados de modelo;
number = número de arquivos;
levtype = nível: sl = surface level e pl = pressure level;
date = intevalo dos dados que você irá baixar;
time = pode ser 00, 06, 12, 18 horas;
origin = este parâmetro se refere ao tigge (qual modelo você irá baixar);
type = dado: an = análise e pf = previsão;
param = variável;
area = para os dados do tigge você pode escolher a área para os dados do ERA-Interim você baixa para o globo todo;
grid = espaçamento de grade;
target = nome do arquivo que será salvo. (Sempre com a extensão .grib o sistema ainda não possibilita baixar em outra extensão, mas esta em fase de construção baixar em Netcdf por exemplo);

Quando eu cheguei neste ponto eu decidi testar este script para baixar os dados de vento zonal para um determinado período e não consegui encontrar o nome exato para o param(variável) foi ai que encontrei esta informação diretamente no site do ECMWF.

Quando você for no site para download do ERA-Interim (aqui) escolha o período, nível (se for o caso), variável e etc exatamente como fazia anteriormente. No entanto, quando você for clicar para fazer o download (Retrieve Grib ou Retrieve Netcdf) você deverá clicar em Views the MARS request e as informações dos seus dados necessárias para você modificar o seu script aparecerão ai é só modificar e sair para o abraço.

Eu fiz um teste para os seguintes dados:

Eu primeiro escolhi no site Pressure Level;
Escolhi janeiro e fevereiro de 2013;
No Select time escolhi (00:00:00, 06:00:00, 12:00:00, 18:00:00)
No Select Level e parameter escolhi U component of wind e level 850
Depois cliquei em Views the MARS request e apareceu o seguinte:

request (desconsidere)

Estimated number of fields: 472 (desconsidere)

retrieve, (desconsidere)
levelist=850,
stream=oper,
levtype=pl,
param=131.128,
dataset=interim,
step=0,
grid=0.75/0.75,
time=00/06/12/18,
date=2013-01-01/to/2013-02-28,
type=an,
class=ei

E o meu script atualizado ficou assim:


#!/usr/bin/env python
from ecmwfapi import ECMWFDataServer


server = ECMWFDataServer()

server.retrieve({
    'dataset' : "interim", (modificado)
    'step'    : "0", (modificado)
    'levelist' : "850", (inserido)
    'stream' : "oper", (inserido)
    'number'  : "all", (removi)
    'levtype' : "pl", (modificado)
    'date'    : "2013-01-01/to/2013-02-28", (modificado)
    'time'    : "00/06/12/18", (modificado)
    'origin'  : "all", (removi)
    'type'    : "an", (modificado)
    'class'    : "ei", (inserido)
    'param'   : "131.128", (modificado)
    'area'    : "70/-130/30/-60", (removi)
    'grid'    : "0.75/0.75", (modificado)
    'target'  : "uwnd.01012013.01022013.grib" (modificado)
    })

No levelist você pode inserir os níveis que precisar. Espero que tenham entendido e principalmente que funcione.

Boa Sorte e Bom Trabalho a Todos.