GrADSスクリプト

GrADSの説明ページ

GrADSの説明は、東北大学の
・http://wind.gp.tohoku.ac.jp/index.php?%B8%F8%B3%AB%BE%F0%CA%F3%2FGrADS%2FGrADS%A4%CETips
・http://wind.gp.tohoku.ac.jp/index.php?%B8%F8%B3%AB%BE%F0%CA%F3%2FGrADS%2FGrADS%A5%B9%A5%AF%A5%EA%A5%D7%A5%C8%A4%CETips
がよくまとまっていておすすめです

まず、GrADSをつかってみる

ここでは、すでにデータセットとして整備されている、観測に基づいた海洋の水温・塩分の気候値World Ocean Atlas 98を書きます。
場所は　http://www.esrl.noaa.gov/psd/data/gridded/data.nodc.woa98.html
水温の月平均の気候値は　ftp://ftp.cdc.noaa.gov/Datasets/nodc.woa98/temperat/monthly/otemp.anal1deg.nc
塩分の月平均の気候値は　ftp://ftp.cdc.noaa.gov/Datasets/nodc.woa98/salinity/monthly/salt.anal1deg.nc
です。

$ wget ftp://ftp.cdc.noaa.gov/Datasets/nodc.woa98/temperat/monthly/otemp.anal1deg.nc

でダウンロードできます。

ターミナル上で

$ grads

としてGrADSを起動します。
Landscape mode? と聞かれますので、何も書かずにリターン(横長のGrADSウインドウ)かnを書いてリターン(縦長のGrADSウインドウ)します。どっちでも大丈夫ですがなにも押さずにリターンして横長の画面にするほうがやりやすいです。

つづいて、先ほどダウンロードした海洋水温ファイルを開きます。.nc(NetCDFファイル)を開くときにはsdfopenで開きます。

ga-> sdfopen otemp.anal1deg.nc

で開きます。

ga-> q ctlinfo

とするとファイルに関する情報が出てきます。

dset otemp.anal1deg.nc
title NODC World Ocean Atlas 1998
undef -9.96921e+36
dtype netcdf
xdef 360 linear 0.5 1
ydef 180 linear -89.5 1
zdef 24 levels 1500 1400 1300 1200 1100 1000 900 800
 700 600 500 400 300 250 200 150 125 100
 75 50 30 20 10 0
tdef 12 linear 00Z01JAN0001 1mo
vars 1
otemp=>otemp  24  t,z,y,x  Ocean Temperature, analyzed mean, 1-deg grid, Monthly
endvars

xdefがx方向の軸(x=1が東経0.5度、xが1大きくなるごとに1度ずつ東へいく)、ydefがy方向の軸(y=1が南緯89.5度、yが1大きくなるごとに1度ずつ北へ行く)、zdefがz方向の軸(1番目が水深1500m、2番目が1400m、最後の24番目が海面0m)、tdefが時間方向の軸(1番目が0001年1月1日0時、しかしこれは気候値なので年には大して意味はない、tが1大きくなるごとに1mo、すなわち1月ずつずれる)、varsが格納されている変数の数でここでは海水温の1個だけ、otempが変数の名前なので、描画する時にはd otempとします。、次の24は鉛直方向の層数、t,z,y,xは定義座標で時間、水深、緯度、経度の4軸で表されることを意味します。
たとえば、以下の図を描画できます。

ga-> d otemp(z=24,t=1) で1月の海面水温の全球分布
ga-> d otemp(z=24,t=8) で8月の海面水温の全球分布
ga-> d otemp(z=1,t=1) で1月の1500m水温の全球分布

ga-> set lat 20 50
ga-> set lon 120 150
ga-> d otemp(z=1,t=1) で日本付近に注目した1月の海面水温分布

ga-> set lat 30
ga-> set lon 140
ga-> set z 1 24
ga-> d otemp(t=8)　で北緯30度東経140度における8月の水温の鉛直分布

ga-> set lat -90 90
ga-> set lon 180
ga-> set z 1 24
ga-> d otemp(t=8)　東経180度における8月の水温の鉛直緯度分布

ga-> set lat -90 90
ga-> set t 8
ga-> set z 1 24
ga-> d ave(otemp,x=1,x=360)　8月の全球東西平均水温の鉛直緯度分布

などがかけます。

画像の保存

ga-> print

とすると画面にかかれているものがgrads.epsに出力されます。ただし起動時にlandscapeモードをnにせずに起動した場合は90度回転された状態ででてくるので

ga-> print -R

とすると画面に出ている向きに出力されます。

ga-> !convert grads.eps grads.png

とすると.pngファイルに変換されます。

アスキーの数字ファイルをGrADSで読む

海洋の気候値PHCをGrADSでかくことにします。PHCとは海洋の観測と再解析に基づく気候値。http://psc.apl.washington.edu/nonwp_projects/PHC/Data3.html。まず、Temperatureのannual fieldをダウンロードしてください。

$ wget http://psc.apl.washington.edu/nonwp_projects/PHC/Data3/Temp00_p3.obj.gz
$ gunzip Temp00_p3.obj.gz 　　　（で解凍する）

Temp00_p3.objをメモ帳などで開くと、数字がひたすら並んでいます。これはfortranの書式10f8.4のformatで経度方向に360個、緯度方向に180個、鉛直方向に33個の水温データが並んでいます。これを、fortranプログラムを使ってGrADSで読むための.grdファイルに変換します。

↓fortranプログラムの例

program bedelev
implicit none

integer :: i,j,k,n,nx,ny,d,nodata,nrec
real(4) :: a(360,180,33)
real(4) :: b(2138400)
real(4) :: depth(33)
character(16) :: c
open(10,file="Temp00_p3.obj")
open(90,file="Temp00.gmt")
open(99,file="Temp00.grd",form="unformatted",access="sequential",convert="big_endian")
open(98,file="Temp00.ctl")
depth(1)=0.0
depth(2)=10.0
depth(3)=20.0
depth(4)=30.0
depth(5)=50.0
depth(6)=75.0
depth(7)=100.0
depth(8)=125.0
depth(9)=150.0
depth(10)=200.0
depth(11)=250.0
depth(12)=300.0
depth(13)=400.0
depth(14)=500.0
depth(15)=600.0
depth(16)=700.0
depth(17)=800.0
depth(18)=900.0
depth(19)=1000.0
depth(20)=1100.0
depth(21)=1200.0
depth(22)=1300.0
depth(23)=1400.0
depth(24)=1500.0
depth(25)=1750.0
depth(26)=2000.0
depth(27)=2500.0
depth(28)=3000.0
depth(29)=3500.0
depth(30)=4000.0
depth(31)=4500.0
depth(32)=5000.0
depth(33)=5500.0



do i=1,2138400,10
read(10,'(10f8.4)') b(i),b(i+1),b(i+2),b(i+3),b(i+4),b(i+5),b(i+6),b(i+7),b(i+8),b(i+9)
write(*,*) i
end do

n=1
do k=1,33
do j=1,180
do i=1,360
a(i,j,k)=b(n)
n=n+1
end do
end do
end do


do k=1,33
do j=1,180
do i=1,360
write(90,'(f9.2,$)') i-0.5
write(90,'(f9.2,$)') j-90.5
write(90,'(f9.2,$)') depth(k)
write(90,'(f9.2)') a(i,j,k)
end do
end do
end do
nrec=1

write(99) a  !書き込み

write(98,*) "DSET ^Temp00.grd"
write(98,*) "TITLE XXX"
write(98,*) "OPTIONS big_endian"
write(98,*) "UNDEF -99.0"
write(98,*) "XDEF 360 linear 0.5 1.0"
write(98,*) "YDEF 180 linear -89.5 1.0"
write(98,*) "ZDEF 33 levels 0 10 20 30 50 75 100 125 150 200 250 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1750 2000 2500 3000 3500 4000 4500 5000 5500"
write(98,*) "TDEF 1 linear 21Z24OCT2013 180mn"
write(98,*) "VARS 1"
write(98,*) "to 33 0 temperature"
write(98,*) "ENDVARS"

stop
end program

を回すとTemp00.grdとTemp00.ctlが生成される。GrADSでTemp00.ctlを開くと読める。(Temp00.gmtというファイルも出てきますがこれはGMT用です、GrADSでは使用しません)

open Temp00.ctl

あとはWorld Ocean Atlasの場合と同様です。

real(4) :: a(360,180,33) ←単精度実数,x方向に360,y方向に180,z方向に33を指定
open(99,file="Temp00.grd",form="unformatted",access="sequential",convert="big_endian") ←form="unformatted"書式なし,access="sequential"
write(99) a　←これでTemp00.grdに書き込む

がポイントです。

GrADSで書いたデータをバイナリに出力する

ファイルをopenしたあと、以下のようなgrads scriptを動かす。

ga-> set x 1 144 #xmaxは自分で指定する
ga-> set y 1 73 #ymaxも自分で指定する
ga-> set z 1 1 #zmaxは自分で指定する
ga-> set t 1 360 #tmaxも自分で指定する
ga-> set gxout fwrite
ga-> set fwrite -be file1.grd #$outfileは自分で指定。little endianなら-beのかわりに-le
d (変数)
disable fwrite

するとfile1.grdが出力される。このfile1.grdは、以下のfortran90を実行することで読むことができる。

integer,parameter :: xmax=144,ymax=73,zmax=1,tmax=360
real(4) :: a(xmax,ymax,zmax,tmax),b(xmax,ymax,zmax,tmax),c(xmax,ymax,zmax,tmax)
integer :: x,y,z,t,n,i,j,k,nrec

open(11,file="file1.grd",form='unformatted',access='direct',recl=4*xmax*ymax)

nrec=0
do t=1,tmax
 do z=1,zmax
  nrec=nrec+1
   read(11,rec=nrec) ((a(x,y,z,t),x=1,xmax),y=1,ymax)
 end do
end do

a(x,y,z,t)に値が入っています。

時間方向にも軸がある4次元データのascii→gradsへの変換

integer,parameter :: nx=360,ny=180,nz=44,nt=12 !←自分で数字をいれる。この値は1でもかまいません。
integer :: nrec
real(4) :: var(nx,ny,nz,nt)
open(51,file="var.grd",form="UNFORMATTED",access="SEQUENTIAL",convert="BIG_ENDIAN")

nrec=1
do t=1,nt
 do k=1,nz
write(51,rec=nrec) ((var(i,j,k,t),i=1,nx),j=1,ny)
nrec=nrec+1
 end do
end do

とします。　最終的なnrecの値は、「時刻の数×z軸」の数になっている必要があります。　

3.GrADS形式の、１つのgrdファイルに複数の変数のデータを書き込む

すぐやり方を忘れるので自分用メモ。fortranで、

open(51,file="file.grd",form="UNFORMATTED",access="DIRECT",recl=4*nx*ny)
nrec=1
do t=1,nt
write(51,rec=nrec) ((d(i,j,t),i=1,nx),j=1,ny)
nrec=nrec+1
write(51,rec=nrec) ((s(i,j,t),i=1,nx),j=1,ny)
nrec=nrec+1
end do

とすればできる。最終的なnrecの値は、「時刻の数×変数の数×z軸」の数になっている必要があります。

GrADSスクリプト

GrADSのコマンドをいちいち手書きで打ち込んでいくのは効率が悪いので、すべての処理を1つのファイルにまとめて実行します。
.gsファイルにコマンドを打ち込みます。

ga-> !emacs script1.gs &

としてメモ帳を開きます。ここにコマンドを打ち込むのですが、以下のようにコマンドの初めと終わりを'(クオーテーション)で囲みます。

'sdfopen otemp.anal1deg.nc'
'set lat 20 50'
'set lon 120 150'
'd otemp(z=1,t=1)'

このようなファイルを作成した後、

ga-> run script1.gs

とすると

ga-> sdfopen otemp.anal1deg.nc
ga-> set lat 20 50
ga-> set lon 120 150
ga-> d otemp(z=1,t=1)

と同じ効果になります。

サンプル：2013年1月のジオポテンシャル高度平年差

データはncep ncar reanalysisから。ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis.dailyavgs/

サンプルスクリプト

'reinit'
'sdfopen hgt.day.1981-2010.ltm.nc'
'set lev 500'
'set map 3 1 5'

*calendar set
start=01JAN2013
starti=1
end=31JAN2013
endi=31

'ave1=ave(hgt.1,t='starti',t='endi')'


'sdfopen hgt.2013.nc'
'set lat 20 90' ; 'set lon 0 360' ; 'set mproj nps'
'set gxout shaded'
'run redblue.gs'
'set clevs -180 -150 -120 -90 -60 -30 0 30 60 90 120 150 180'
'set ccols 20 21   22   24   26  28 29 31 32 34 36 38 39 40' 
'd ave(hgt.2,time=00Z'start',time=00Z'end')-ave1'
'cbarn'
'set gxout contour'
'set cstyle 1'
'set ccolor 1'
'set cmin 4860'
'set cmax 6000'
'set cint 60'
'd ave(hgt.2,time=00Z'start',time=00Z'end')'
'draw title 'start'-'end' 500hPa Geopotential Height'

GrADSがうまくいかないぞ

値がおかしい

1e+34みたいな値になっていたらendianの問題である可能性が高い。endianをbigならlittleにかえたりその逆を試す。またはexport F_UFMTENDIAN=bigなどを試す。あとはファイル書き込み時にreal(8)にしていないか。

変数を書くとエラーになる

.ctlファイルで、VARSの次の行、変数の次の数字は鉛直層数になっていなければいけない。
自分で.grdファイルを作った場合なら、配列の配置が間違っている

そもそもファイルが開けない

.ncファイルなら必ずsdfopenできるわけではない。COARSE規約でないNetCDFはsdfopenできません。この場合、自分で.ctlファイルを使って.ncファイルを開きます。NetCDFをGrADSで開く項目を見てください。
自分で.grdファイルを作った場合なら、.ctlファイルに誤りがある。

.ncファイルの情報が分からない（変数単位が分からないとか）

ncdump -c (.ncファイル) をtmpなどにリダイレクトしてください。この中に情報が入っています。NetCDFをGrADSで開く項目を見てください。

やってみたとおりにしてもうまくいかないぞ

計算機環境に依存するところは出てきてしまうので試行錯誤するしかないこともあります。

back