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!

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!!

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!!

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!!

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!

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!

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!