公 開 されているデータのページ まずは sea level pressure の ltm( 長 期 間 平 均 )データをダウンロードしてみよう リンクにマウスを 合 わ

GrADSで図を描く練習

GrADS講習会 第一回 09 年4 月作成 14年1月修正 Y. Kamae http://seesaawiki.jp/ykamae_grads-note/ http://researchmap.jp/muuokv4p4-1772021/#_1772021 2

GrADSに関して詳しいページ

http://seesaawiki.jp/ykamae_grads-note/ GrADS-Note GrADSコマンドやスクリプトのメモ、インストールについて http://wind.geophys.tohoku.ac.jp/index.php?%B8%F8 %B3%AB%BE%F0%CA%F3%2FGrADS 東北大学大学院理学研究科 流体地球物理学講座 公開情報/GrADS 必ずブックマークしておきたい http://radar.sci.hokudai.ac.jp/~yoshida/study/GrADS/Command/ GrADS コマンドリファレンス 吉田さんによるコマンドリファレンス日本語訳 http://hydro.iis.u-tokyo.ac.jp/~kei/?IT%20memo%2FGrADS%20memo IT memo 芳村さんによるやや上級者向け情報 http://mausam.hyarc.nagoya-u.ac.jp/%7Ehatsuki/grads/gradsmanu.html GrADS リファレンスマニュアル 藤波さん作成マニュアル http://www.iges.org/grads/gadoc/gadocindex.html GrADS Documentation Index

GrADS本家のコマンド説明（英語） 3

1. NCEP/NCAR 再解析データ

2. wget

データのダウンロード

3. sdfopen

NetCDFファイルを開く

4. q file

ファイルの情報を調べる

5. d 変数名 図を描く

6. xanim

アニメーションを描く

今日の内容 4 研究を行なうには、データが必要。 研究に使うデータはどこにある？ ・観測 自分で収集する ・解析 誰かが作ったデータを使う（陸域センター、気象庁、研究所、・・・） アメリカの大気研究センター、環境予測センターが管理している全球のデータセット 1949年～現在まで ○世界で一番降水が多い地域は？ ○日本の平年の冬の気圧配置 ○過去50年の気温の変化傾向 いろんなことを調べることができます 今回は、各月の平年の値（気候値）のデータについて見てみましょう

NCEP/NCAR

再解析データ

NCEP/NCAP 再解析データ

5 http://www.cdc.noaa.gov/cdc/data.ncep.reanalysis.derived.surface.html

公開されているデータのページ

まずは sea level pressure の ltm（長期間 平均）データをダウンロードしてみよう リンクにマウスを合わせ、右クリック* →リンクのURLをコピー →ターミナルを開き、マウスの中クリック* URLがペーストされる 入力したURLの前に、wget と入力 つまり

wget

データをダウンロード

する

wget ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis.derived/surface/slp.mon.ltm.nc を、実行 ファイルがダウンロードされる ls と打って、ファイルがダウンロードされたか確認

wget ファイルの場所

ファイルのダウンロード *2ボタンマウスでは左右同時クリック、 Macではoption+クリック *Macではcontrol+クリック 6

sdfopen

nc

データを

GrADS

で開く

と打ち、GrADSを起動させる nc (NetCDF) 形式のファイルは、 sdfopen ファイル名 で開くことができます つまり grads sdfopen slp.mon.ltm.nc と入力 ※プログラミングなどでよく使うバイナリ（bin）形式のデータを開くときは、 コントロール（ctl）ファイルを作る必要があります。 名前 拡張子 GrADSで開くコマンド

NetCDF

.nc

sdfopen ~.nc

コントロールファイル

_.ctl

_{open ~.ctl}

7

q file

データの情報を

調べる

q file と打つと、今開いたデータの情報 を調べることができる

File 1 : Longterm Monthly Mean SLP NCEP Reanalysis Descriptor: slp.mon.ltm.nc

Binary: slp.mon.ltm.nc Type = Gridded

Xsize = 144 Ysize = 73 Zsize = 1 Tsize = 12 Number of Variables = 1

slp0 -999 Monthly Longterm Mean of Sea Level Pressure

1番のファイルとして認識 ファイルについての説明 ファイルの名前 東西、南北、高さ、時間方向の データの数 変数名 d slp d のあとに変数名を入力すると、図が描ける 最初に入っていたデータ（1月）の図 つまり 1968 - 1996 年平均の1月の海面気圧が描かれる

d

変数名

図を描く

8

xanim

アニメーション

set t 2 d slp 時間を2ステップ目（2月）に設定 図を描く c 一度図を消して set t 1 12 1月から12月まで d slp アニメーションで見る ←早すぎてよくわからない・・・ c set t 1 12 xanim ‐pause slp 1月から12月までのアニメーション。 図をクリックするたびにひと月進む。 上に時間も表示されてわかりやすい！ －－練習－－ slp以外の好きなファイルをダウンロードして、図を描いてみよう！ ※降水について描いてみたい場合は ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis.derived/surface_gause/prate.mon.ltm.nc をwget http://seesaawiki.jp/ykamae_grads-note/ の GrADSスクリプトライブラリ 参照 ※ ※

エルニーニョを描く

GrADSとGrADSスクリプト

GrADS講習会 第二回 09 年4 月作成 14年1月修正 Y. Kamae http://seesaawiki.jp/ykamae_grads-note/ 2 局地的な観測データはテキスト形式、衛星データはhdf形式で公開されているものが多いが、 全球データはNetCDF形式のものが多い →GrADSが大活躍 今回の使用データ NOAA（アメリカ海洋大気圏局 ）サイト http://www.esrl.noaa.gov/psd/data/gridded/data.noaa.ersst.html にある、ERSST 公開データ sst.mnmean.nc

9798elnino.gs cor_reg.sh fwrite-sst.gs mon-ave.gs mon-fig.gs nino3ano.gs sst1980-1999.ctl sst_mul.gs sst.mnmean.nc を GrADSで開いてみよう cd 02grads unzip 02grads.zip wget ftp://ftp.cdc.noaa.gov/Datasets/noaa.ersst/sst.mnmean.nc にある ”02grads.zip” をダウンロード、作業ディレクトリに置く http://researchmap.jp/muwnuoeel-1772021/#_1772021 xtermを起動（MacならX11、cygwinならstartxwin.sh etc） 3 grads

File 1 : NOAA Extended Reconstructed SST V2 Descriptor: sst.mnmean.nc

Binary: sst.mnmean.nc Type = Gridded

Xsize = 180 Ysize = 89 Zsize = 1 Tsize = 1850 Number of Variables = 1

sst0 -999 Monthly Means of Sea Surface Temperature 1番目に開いたファイルなので、file 1 x, y, z, t方向にどれだけデータがあるか 変数名 GrADSでデータを開いたときは、まずq fileをして、ファイルの情報を確認しよう ここで変数名が何かわかるので、dで図を描く ファイルの情報を見る NetCDFファイルを開く GrADSを起動、Enter2回、白画面クリック q file sdfopen sst.mnmean.nc 変数名sstについての図を描く d sst 4 GrADSに関しては、学群実験のテキストや 東北大学大学院理学研究科 流体地球物理学講座 公開情報/GrADS http://wind.geophys.tohoku.ac.jp/index.php?%B8%F8%B3%AB%BE%F0%CA%F3%2FGrADS を参照すれば、大抵のことはできるようになります などなど 初期状態に戻す カラーバーを表示 新しい設定で図を描き 背景を白に変更 シェードを描こう 図を消す reinit cbar d sst c

set display color white set gxout shaded c

5 領域平均： その変数について、東西A~B°、南北C~D°の領域で平均した値

領域平均

北緯4度から南緯4度，西経150度から90度の海域はNino3と呼ばれ，気象庁は Nino3の月平均海面水温を用いてエル・ニーニョを定義しています。 この領域で平均した値を使って、エルニーニョの強弱の時系列を書いてみよう

aave(変数名, lon=A,lon=B,lat=C,lat=D)

sdfopen ～.nc 0 360 0 -90 90 210 270 ↑領域平均の値をaと定義し、aを描いた t=1 で 25.4629 となる d a a=aave(sst,lon=210,lon=270,lat=-4,lat=4) 6 reinit

時系列

水平方向に幅のあるデータは、何も設定しないと水平分布図が描かれます。 このとき水平方向の幅をなくし、代わりに時間方向に幅を持たせると、時系列図 が描かれます。 ←水平方向の幅をなくす 10年間のNino3 SSTの時系列が描かれた ←時間方向に幅を持たせる d a a=aave(sst,lon=210,lon=270,lat=-4,lat=4) set t 1 120 set y 1 set x 1 sdfopen ～.nc 7

時間平均、気候値、偏差

d ave(変数名,t=a,t=b) で、時間a~bの期間の平均値が描かれるのはおなじみ 1980年1月から、1999年12月までの、偏差時系列を描いてみよう reinit 82/83年、97/98年が強いエルニーニョ年だったことがわかる ←1980年1月～1999年12月 ←偏差＝生の値－気候値（20年間平均値） d a-clim a=aave(sst,lon=210,lon=270,lat=-4,lat=4) set t 1513 1752 set y 1 set x 1 sdfopen ～.nc set t 1 12 clim=ave(a,t+1512,t=1752,12) modify clim seasonal set t 1513 1752 ←1月から12月までの ←各月ごとの20年間平均値（気候値）を計算 ←気候値climを12か月周期のデータにする 8

GrADSスクリプト

今までやってきたことを、レシピのように書き留めて、一発で図を描かせることができる nino3ano.gs 偏差時系列を図化 'reinit' 'sdfopen sst.mnmean.nc' 'set x 1' 'set y 1' 'set t 1513 1752' 'a=aave(sst,lon=210,lon=270,lat=-4,lat=4)' 'set t 1 12' 'clim=ave(a,t+1512,t=1752,12)' 'modify clim seasonal' 'set t 1513 1752' 'd a-clim' 'print nino3ano.eps' 9798elnino.gs 強Elnino時の偏差を図化 'reinit' 'sdfopen sst.mnmean.nc' 'set t 1 12' 'clim=ave(sst,t+1512,t=1752,12)' 'modify clim seasonal' 'set t 1' 'd ave(sst-clim,t=1727,t=1729)' 'print 9798elnino.eps' で図ができているか見てみよう run ～.gs で実行 gradsを起動、 gs ～.eps gradsを終了したのち、

ctlファイルとGrADSスクリプト

binファイルへの書き出し、

ループで作図

GrADS講習会 第三回 09 年4 月作成 14年1月修正 Y. Kamae http://seesaawiki.jp/ykamae_grads-note/ 2 第二回で扱ったデータを使います ・02gradsディレクトリがある人は、そこへ移動 ・ ない人は、第二回 2ページ目のコマンドを実行

nc

gs

png, eps他

データ切り出し 計算など gs GrADSスクリプト nc NetCDFファイル bin バイナリデータ ctl コントロールファイル png, eps 画像の形式

bin

bin

データ 格子情報

ctl

作図 データと格子情報を 含む元データ binデータの格子情報 をGrADSに伝える GrADSに作業をさせる スクリプト

ファイルの種類と

役割

3 sst.mnmean.nc fwrite-sst.gs sst1980-1999.bin _{sst1980-1999.ctl} mon-ave.gs mon-ave.bin _mon-ave.ctl mon-fig.gs mon-1.png ・・・ mon-12.png 154年×12ヶ月の データ・格子情報 ①20年間のデータだけ切り出す ②各月の気候値（平年値）を作成 ③各月の図を描く

作業の流れ

① ② ③ 20年×12ヶ月 データ 気候値×12ヶ月 格子情報 スクリプト 4

binファイルへの書き出し、ctlファイル

fwrite-sst.gs 必要な20年間だけ取り出してbinに変換 ‘reinit’ ‘sdfopen sst.mnmean.nc’ ‘set x 1 180’ ‘set y 1 89’ ‘set t 1513 1752’

‘set fwrite -le sst1980-1999.bin’ ‘set gxout fwrite’

‘d sst’ ‘disable fwrite’

sst1980-1999.ctl

binをGrADSで読むためのファイル dset ^sst1980-1999.bin

titile monthly sst 1980-1999 NOAA ERSST undef 32767

options 365_day_calendar little_endian xdef 180 linear 0 2.0

ydef 89 linear -88.0 2.0 tdef 240 linear 1jan1980 1mo zdef 1 levels 1

vars 1 sst 1 0 x,y sst endvars

gradsを起動させ、run fwrite-sst.gs で実行

reinit し、open sst1980-1999.ctl で開き、 d sst で描かれるのは1980年1月のデータ データをファイルに 書き出す データ元サイト やq fileの情報 をもとに記述 読み込みたい binファイル ctlファイルの書式については、 http://seesaawiki.jp/ykamae_grads-note/ →”GrADS ctlファイルの作成” などを参考に

5 mon-ave.gs 20年（240カ月）のデータを使って、各月の気候値（平年値）を作る 'reinit' 'open sst1980-1999.ctl' 'set x 1 180' 'set y 1 89' 'set gxout fwrite'

'set fwrite -le mon-ave.bin' t=1 while(t<=12) 'set t 't'' 'd ave(sst,t+0,t=240,12)' t=t+1 endwhile 'disable fwrite'

ループを使って書き出し

基本的にfwrite-sst.gs と同じ クオーテーションで囲ま ない場所は、grads自体 には命令されない t の値が1ずつ増えて いく ＝tのループが1～12 の範囲で回る※

grads -lc mon-ave.gs で実行

grads -lc ～.gs

～.gsの実行 ※この場合はset t 1 12でも同じだが、後述の図の書き出しなどの際はループを使う必要がある 6 cp sst1980-1999.ctl mon-ave.ctl vi mon-ave.ctl できたbinファイル「mon-ave.bin」を読むctlファイルを作ろう（以前のctlファイルをコピーして書き換 えよう） mon-ave.ctl dset ^mon-ave.bin

title monthly average sst 1980-1999 NOAA ERSST undef 32767

options 365_day_calendar little_endian xdef 180 linear 0 2.0

ydef 89 linear -88.0 2.0 tdef 12 linear 1jan1980 1mo zdef 1 levels 1 vars 1 sst 1 0 x,y sst endvars 赤の部分を書き換える だけでよい titleは書き変えなくても 問題ない ・ 対応するbinファイル名 ・ 変更した空間、時間の情報 を書き換えればよい 今回の場合、読みこむbinファイルがmon-ave.binになったこと、時間の総数が240から 12になったこと、の2点だけ変更すれば、あとはGrADSを起動させて open mon-ave.ctl d sst で、図にできることが確認できる 7 'reinit' 'open mon-ave.ctl' 'set lon 40 240' 'set lat -40 60' 'set display color white' t=1

while(t<=12) 'c' 'set t 't'' 'set gxout shaded' 'set rgb 20 0 0 255' 'set rgb 21 100 100 255' 'set rgb 22 200 200 255' 'set rgb 23 255 255 255' 'set rgb 24 255 200 200' 'set rgb 25 255 100 100' 'set rgb 26 255 0 0' 'set clevs 5 10 15 20 25 30' 'set ccols 20 21 22 23 24 25 26' 'set grads off'

'd sst' 'cbar'

'set gxout contour' 'set clopts 1 2 0.18'

'set clevs 0 2.5 5 7.5 10 12.5 15 17.5 20 22.5 25 27.5 30 32.5'

'set cthick 4' 'set ccolor 1' 'set grads off' 'd sst'

'draw title NOAA ERSST clim sst <1980;1999> mon='t'' 'printim mon-'t'.png' t=t+1 endwhile mon-fig.gs まずは、grads -lc mon-fig.gs で実行 「各月のアジア域の図が、低温青、高温赤のシェードと コンターで描かれる」 quit で終了し、ls してみると、12個の図ができている ことがわかる アカウント名ディレクトリ→02grads に移動してみると、それぞれの図があるはず それぞれの行の役割を、東北大のGrADSページを見て調べ てみよう 8 「気候値SSTの季節変化」の図の完成！

カラーパレットの設定

図の書き出し

GrADS講習会 第四回 09 年4 月作成 14年1月修正 Y. Kamae http://seesaawiki.jp/ykamae_grads-note/ 2 第二回・第三回で扱ったデータを使います ・02gradsディレクトリがある人は、そこへ移動 ・ない人は、第二回 2ページ目のコマンドを実行 3 grads GrADSを起動、Enter2回、白画面クリック

sdfopen sst.mnmean.nc NetCDFファイルを開く

q file ファイルの情報を見る（一番左下に変数名が表示される）

d sst 変数名sstについての図を描く

c 図を消す

set gxout shaded シェードを描く

set display color white 背景を白に変更

c d sst 新しい設定で図を描く cbar カラーバーを表示 c 図を消す

カラーパレット（初期設定）

コンター、シェードを描くとき、カラーパレットの初期設定はレインボーになっている 4

カラーパレットの設定

set rgb 20 200 200 200 set rgb 21 100 100 100 set rgb 22 0 0 0 set clevs 15 25 set ccols 20 21 22 d sst cbar この7行を打ち込むと、↓の3色の図になる

set rgb 20 200 200 200

色番号※_{＝20を赤＝200、緑＝200、青＝} 200に設定

set clevs 15 25

境界の値を「15」「25」に設定

set ccols 20 21 22

clevsで設定した境界の値に従って、色番 号をあてはめる ※色番号 1~15番は既に設定されている 自分で色を設定したい場合は、16以降を使う 値 ～15 15～25 25～ 番号 20番 21番 22番 赤緑 青 200 200 200 100 100 100 0 0 0

5

rgbの設定

光の3原色を使って色を設定する

r

g

b

黒

0

0

0

赤

255 0

0

緑

0

255 0

青

0

0

255

黄

255 255 0

紫

255 0

255

青緑

0

255 255

白

255 255 255

R

G

B

http://ja.wikipedia.org/wiki/RGB

その他の色の例は http://www.scollabo.com/banban/lectur/websafe.html http://seesaawiki.jp/ykamae_grads-note/ → GrADSスクリプトライブラリ 参照 6

colorコマンド

colorというコマンド※_{を使えば、いちいち設定} しなくても自動で色をつけてくれる c color 0 30 2.5 d sst cbar 0から30まで、2.5きざみで、低いほど寒色系、 高いほど暖色系のカラーパレットを作成

color min max int

※正確にはコマンドではなくスクリプト。color.gsという GrADSスクリプトを指定の場所に置いておくことで、 「color」と打ちこんだときにそのスクリプトが起動される。具 体的にはgrads本体があるディレクトリの一つ上のlib/の中 など minからmaxまで、int間隔で設定 （コンターを引くときも使える） 7

色を塗るときの注意

IPCC 2007 Ch.11より 降水量の 将来変化 地上気温の 将来変化 夏季平均 温かい、熱い イメージ 冷たい、寒い イメージ 気温の変化 湿潤な イメージ 乾燥した イメージ 降水量の変化 寒色系 暖色系 茶色、赤 緑、青

color min max int -kind blue->white ->yellow->red->darkred

color min max int -kind saddlebrown->white->green

描く物理量によってカラーパレットを変えると、イメージが伝わりやすい 8

図の書き出し

画面に表示させた図をファイルにするには？ png形式にしたい場合 printim ～.png eps形式にしたい場合 print ～.eps ※Ver.1.8以前のGrADSでは enable print hoge.gx print

disable print

!gxeps –c –i hoge.gx –o hoge.eps

png形式の図を表示（ImageMagick） display ～.png eps形式の図を表示（ghostscript） gs ～.eps ※Ver.1.9以上のGrADSのみ 作成した図を表示するには？ pptに載せる程度であれば、pngで十分 epsに対応した環境であれば、epsを使用。 png（や、jpg、gif）を使いたい場合は、まずepsで図を出力し、 convertで高画質変換する。 GrADS内で print hoge.eps コマンドラインで

convert -density 600 -resize 600 -rotate 90 +antialias hoge.eps hoge.png

-densityで解像度を指定、-resizeでサイズを指定、 -rotateで 90度回転させ、 +antialiasで邪魔な横線が入らないように pngに変換 display hoge.png で確認。printimで作成したときより、画質が高いことがわかる。 図の解像度を上げたい

GrADS作図オプション

見やすい図に仕上げる

GrADS講習会 第五回 09 年4 月作成 14年1月修正 Y. Kamae http://seesaawiki.jp/ykamae_grads-note/ 2

タイトル表示draw title CMAP Pr & NCEP/NCAR 850hPa Wind JJA xy軸ラベル間隔 set xlint 10 set ylint 10 xy軸ラベル文字 色（初期1）、太さ4、大きさ0.12 set xlopts 1 6 0.19 set ylopts 1 6 0.19 コンターラベル文字 色（初期1）、太さ1、大きさ0.09 set clopts 1 5 0.22 ラベルの表示・非表示 set clab off set clab on

1. 文字

http://radar.sci.hokudai.ac.jp/%7Eyoshida/study/GrADS/Command/axis.html その他：GrADSコマンドリファレンス 軸ラベル・アノテーション 3 コンターを引く値 set clevs 5 15 20 コンターの太さ set cthick 5 コンターの色 set ccolor 1 ラベルの非表示 set clab off シェードを等間隔表示

最小値 最大値 間隔 色のグラデーション color 2.5 15 2.5 –kind white->green カラーバーの表示 大きさ 縦0/横1 x方向中心位置 y方向中心位置 cbarn 1 0 5.5 0.18

2. コンター・シェード

http://radar.sci.hokudai.ac.jp/%7Eyoshida/study/GrADS/Command/graph_opt.html その他：GrADSコマンドリファレンス グラフィックオプション コンター：set gxout contour シェード：set gxout shaded

色の設定、間隔の設定については第四回参照 4 ベクトルを数点間隔で描く d skip(umask,2);skip(vmask,2) 弱いベクトルの非表示 uとvの√二乗和が～以下のとき非表示 umask=maskout(u,mag(u,v)-2) vmask=maskout(v,mag(u,v)-2)

3. ベクトル

ベクトル：set gxout vector

d 東西風;南北風

http://radar.sci.hokudai.ac.jp/%7Eyoshida/study/GrADS/Command/graph_opt.html

5 'reinit'

'sdfopen precip.mon.ltm.nc' 'set gxout shaded' 'set display color white' 'c' 'set vpage 0.1 11.0 0.3 8.5' 'set lon 40 160' 'set lat -20 50' 'set xlint 10' 'set ylint 10' 'set xlopts 1 6 0.19' 'set ylopts 1 6 0.19' 'set grads off'

'color 2.5 15 2.5 -kind white->green' 'd ave(precip,t=6,t=8)' 'cbarn 1 0 5.5 0.18' 'close 1' 'sdfopen uwnd.mon.ltm.nc' 'sdfopen vwnd.mon.ltm.nc' 'set vpage 0.1 11.0 0.3 8.5' 'set lon 40 160' 'set lat -20 50' 'set z 3' 'set xlint 10' 'set ylint 10' 'set xlopts 1 6 0.19' 'set ylopts 1 6 0.19' 'set grads off'

'uave=ave(uwnd.1,t=6,t=8)' 'vave=ave(vwnd.2,t=6,t=8)' 'umask=maskout(uave,mag(uave,vav e)-2)' 'vmask=maskout(vave,mag(uave,vav e)-2)'

'set mpdraw off'

'd skip(umask,2);skip(vmask,2)' 'close 2' 'close 1' 'sdfopen precip.mon.ltm.nc' 'set vpage 0.1 11.0 0.3 8.5' 'set lon 40 160' 'set lat -20 50' 'set xlint 10' 'set ylint 10' 'set xlopts 1 6 0.19' 'set ylopts 1 6 0.19' 'set gxout contour‘ 'set clevs 5 15 20' 'set cthick 5' 'set ccolor 1' 'set clab off' 'set grads off' 'set mpdraw off' 'd ave(precip,t=6,t=8)' 'set clopts 1 5 0.22' 'set clevs 10' 'set cthick 5' 'set ccolor 1' 'set clab on' 'set grads off' 'set mpdraw off' 'd ave(precip,t=6,t=8)'

'draw title CMAP Pr & NCEP/NCAR 850hPa Wind JJA'

'print sakuzu.eps' 'printim sakuzu.png' 格子情報が異なるため、降水と風は同時に開かず、閉 じてからsetをし直し描く 背景の地図はdのたび色が変わってしまうので、2回目 以降は上書きしない ①降水のシェードを描く ②風のベクトルを描く ③降水のコンター（ラベルなし→ラベル あり）を描く 6

4. 折れ線

折れ線の太さ、スタイル、色 set cthick 9 set cstyle 2 set ccolor 8 数値軸の範囲の設定 set vrange 1 10 気候値の図などで、年を消した いとき

set tlsupp year マーカーの大きさ、スタイル set digsiz 0.15 set cmark 4 http://radar.sci.hokudai.ac.jp/%7Eyoshida/study/GrADS/Command/graph_opt.html その他：GrADSコマンドリファレンス グラフィックオプション http://radar.sci.hokudai.ac.jp/%7Eyoshida/study/GrADS/Command/color.html 7 'reinit' 'sdfopen precip.mon.ltm.nc' 'set x 1' 'set t 1 12'

'set display color white' 'set gxout shaded' 'c' 'set vpage 0.1 11.0 0.1 8.5' 'set vrange 1 10' 'set xlopts 1 6 0.19' 'set ylopts 1 6 0.19' 'set digsiz 0.15' 'set cthick 9' 'set cstyle 2' 'set cmark 4' 'set ccolor 8' 'set tlsupp year' 'set grads off'

'd tloop(aave(precip,lon=70,lon=100,lat=0,lat=20))' 'set cmark 2'

'set cstyle 1' 'set ccolor 9' 'set tlsupp year' 'set grads off'

'd tloop(aave(precip,lon=120,lon=150,lat=-20,lat=0))' 次ページにつづく

4. の作図スクリプト

南アジア 70~100E EQ~20N で定義 北オーストラリア 120~150E 20S~EQ ↑6~8平均－12~2平均の降水量分布 南北半球の乾季と雨季の交替 8 ※使用データ NCEP/NCAR Reanalysis 長期間月平均データ 東西風 ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis.derived/pressure/uwnd.mon.ltm.nc 南北風 ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis.derived/pressure/vwnd.mon.ltm.nc

CMAP(CPC Merged Analysis of Precipitation) standard長期間月平均データ 降水量 ftp://ftp.cdc.noaa.gov/Datasets/cmap/std/precip.mon.ltm.nc http://www.cdc.noaa.gov/PublicData/より つづき （凡例の表示） ‘set line 0’ ‘draw recf 8.4 6.5 10.2 7.6’ ‘set line 1 1 6’ ‘draw rec 8.4 6.5 10.2 7.6’

‘cbar_line_set –x 8.7 –y 7.3 –r 0.08 –n 9 –l 1 –c 9 –m 2 –z 0.16 –sn 6 –sz 0.20 –t “Aust”’ ‘cbar_line_set –x 8.7 –y 6.8 –r 0.08 –n 9 –l 2 –c 9 –m 4 –z 0.16 –sn 6 –sz 0.20 –t “Asia”’ 'draw title CMAP Pr South Asia vs North Australia'

'print hoge.eps'

「color」「cbarn」「cbar_line_set」の使用法について

GrADSのウィンドウと

描画の位置

GrADS講習会 第六回 09 年4 月作成 14年1月修正 Y. Kamae http://seesaawiki.jp/ykamae_grads-note/ 2 第二回・第三回で扱ったデータを使います ・02gradsディレクトリがある人は、そこへ移動 ・ ない人は、第二回 2ページ目のコマンドを実行 3 GrADSの起動時、画面の大きさを変えられる

grads -l

grads -p

grads –lg 250x200

大きさ：500 x 400

格子：11.0 x 8.5

大きさ：400 x 500

格子：8.5 x 11.0

大きさ：250 x 200 格子：11.0 x 8.5 landscapeモード ＝横長 portraitモード ＝横長 gradsの起動オプション -l 横長 -p 縦長 -g 大きさ -b バッチモード（黒画面開かない） -c コマンド名 起動後、このコマンドを実行 -help 起動オプション表示

11.0格子

8.5

格子

11

.0

格子

8.5格子

例. GrADSでhoge.gsをrunさせる → grads –lc hoge.gs 4 sst.mnmean.ncを開き

set display color white c set xlopts 1 4 0.22 set ylopts 1 4 0.22 d 変数 x軸とy軸のラベル を大きく ラベルが左端をはみ出してしまう そこで、図の四角枠の範囲を狭め、画面に収まる ようにする c set parea 1.0 10.0 0.0 8.5 d 変数 parea 1.0 10.0 0.0 11.0 0.0 8.5 図の枠の範囲を 設定 図の枠を狭めたので、ラベルも画面に収まる 図の枠外の範囲を決めたいときは、set vpageを使う

5 一度に複数の図を並べることもできる c d 変数(t=1) 画面の下に mul.gsを使うと、もっと楽！ c mul 1 2 1 1 横1x縦2に分割し (1,1)に set parea 0.0 11.0 0.5 4.0 set parea 0.0 11.0 4.5 8.0 d 変数(t=2) t=1 を描き 画面の上に t=2 を描く 0.5 4.0 4.5 8.0 http://seesaawiki.jp/ykamae_grads-note//→ GrADSスクリプトライブラリ mul.gs について d 変数(t=1) mul 1 2 1 2 d 変数(t=2) (1,2)に t=1 を描き t=2 を描く

1

1

2

ただし色や書式の設定を直さないといけないので、面倒 6 'reinit'

'set display color white' 'c'

'set gxout shaded' 'sdfopen sst.mnmean.nc' 'set parea 0.5 8.0 1.4 5.5' 'set xlopts 1 6 0.19' 'set ylopts 1 6 0.19' 'set ylint 30' 'set x 1 180' 'set y 1 89' 'set time aug2000' 'color -5 30 5' 'set grads off' 'd sst'

'cbarn 1.1 0 4.25 0.4' 'set parea 8.5 10.8 1.4 5.5' 'set gxout line'

'a=ave(sst,x=1,x=180)' 'set x 1'

'set xyrev on' 'set ccolor 1' 'set cthick 10' 'set cmark 0'

'set vrange -10.0 30.0' 'set xlint 10'

'set grads off' 'd a'

'print hoge.eps'

'!convert -rotate 90 -density 600 -resize 600 +antialias hoge.eps sst_mul.png' '!rm hoge.eps' 2000年8月のSST分布と東西平均SST分布 画面分割の例 sst_mul.gs 7 8 GrADSの新しいバージョンについての情報 http://seesaawiki.jp/ykamae_grads-note/d/GrADS2%a4%de%a4%c8%a4%e1 GrADS 2.1以降 gxprint 任意のフォーマットの画像ファイルを出力 set rgbで透過色を作成 set tileでドットや格子状のタイル模様を描く Cairoフォントが使用可能 マークが追加 set gxout shade2に対応 グリッド内挿（解像度の変更） lterp

365 day calendarのNetCDFファイルに対応していない z軸を降順に並びかえ

shapefileの読み書き sdfwrite NetCDFの作成

set clab masked コンターラベルを重ならないように描く 画面アスペクト比を変更(V2.0.a9以降) GRIB2形式のデータを読む wgrib2の入手とコンパイル g2ctlの入手 実行 GrADS2 システム要件 GrADS2のインストール など

ASCIIを変換して

GrADSで読む

Fortran + GrADS

09 年4 月作成 14年1月修正 Y. Kamae http://seesaawiki.jp/ykamae_grads-note/ 2 ascii（人間が読めるデータ）をbinary（機械しか読めないデータ）に変換するプログラム program main implicit none integer,parameter:: nn=5 real:: a(nn) integer:: in open(10,file='data1.txt', & access='sequential',form='formatted') open(20,file='data1.bin', & access='direct',form='unformatted',recl=4*nn) read(10,*) (a(in),in=1,nn) write(20,rec=1) (a(in),in=1,nn) end 変数の数を設定 暗黙の型宣言を解除 変数の型を宣言 -0.1 2.5 4.0 1.1 -10.0 という中身のdata1.txtを作り、これをbinaryに変換してみよう 以上のプログラムをtxt-bin.f90という名前で作成し、 g95 txt-bin.f90 ./a.out

キホン

asciiファイルのopen binaryファイルのopen コンパイラはifort 等でもよい 3 http://islscp2.sesda.com/ISLSCP2_1/data/snow_sea-ice_oceans/snow_cover_xdeg/

実践

例としてISLSCP2_1

のNorthern Hemisphere Monthly Snow Cover Extentデータをbinに変換し、 GrADSで読んでみよう！ cd 02grads ls snow_cover_data_1d_1986.csv というデータが1995年まである。_{csv形式（ASCII）なのでlessで読める。} less snow_cover_data_1d_1986.csv LONG,LAT,JAN,FEB,MAR,APR,MAY,JUN,JUL,AUG,SEP,OCT,NOV,DEC,LAND_MASK -38.50,83.50,100,100,100,100,100,100,100,100,100,100,100,100,1 -37.50,83.50,100,100,100,100,100,100,100,100,100,100,100,100,0 ・・・ というような並びになっている。 このデータをfortranで読み、緯度経度時間並びのbinaryに変換してみよう。 4 program main integer:: c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,d real:: a,b

open(10, file='snow_cover_data_1d_1986.csv', &

status='unknown',access='sequential',form='formatted') read(10,'(66a)') read(10,*)& a,b,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,d write(6,*)& a,b,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,d end ASCIIデータをopen ヘッダ行66文字を読む 書式識別子を* にしておくと , などを適当に読み飛ばしてくれる 最初の2つが経度緯度、12個が各月データ、最後が海陸 まずは簡単にf90の基本の書式を書き、データをカンマ区切りで読んでみる。 以上のプログラムをread-csv.f90という名前で作成し、 g95 read-csv.f90 ./a.out -38.50000 83.50000 100 100 100 100 100 100 100 100 100 100 100 100 1 すると と出力される。無事、実数・整数など読み込みに成功したようだ。 コンパイラはifort 等でもよい

5

あとは、宣言文を正確に記述し、最終的に(経度,緯度,時間)の変数としてbinに出力させる。

program main

implicit none !暗黙の型宣言の無効

integer, parameter:: nl=8721,nx=360,ny=60,nmon=12 !整数の定数 real, parameter:: undef=1e+20 !実数の定数

integer:: c(nl,nmon),d(nl) !整数の変数 integer:: iy,il,mon,lat,lon !ループに使う整数 real:: a(nl),b(nl),ilon,ilat !緯度経度(実数) real:: snow(nx,ny,nmon) !書き込むデータ(実数) open(10, file='snow_cover_data_1d_1986.csv', &

status='unknown',access='sequential',form='formatted') open(20, file='test.bin', &

status='unknown',access='direct',form='unformatted',recl=4*nx*ny) read(10,'(66a)') do il=1,nl read(10,*) & a(il),b(il),c(il,1),c(il,2),c(il,3),c(il,4),c(il,5),c(il,6),& c(il,7),c(il,8),c(il,9),c(il,10),c(il,11),c(il,12),d(il) !データcはilとmonの配列にする enddo !il do lat=1,ny do lon=1,nx ilat=lat+23.50 ilon=lon-180.50 do il=1,nl

if((b(il)==ilat).and.(a(il)==ilon)) then !il行目のlatlonがilatilonと一致したとき do mon=1,nmon snow(lon,lat,mon)=real(c(il,mon)) !lonlatmonの関数snow=ilmonの関数c の実数型 enddo !mon goto 100 !一致したときは100へ、一致しなかったときは下へ endif !a,b enddo !il do mon=1,nmon snow(lon,lat,mon)=undef

enddo !mon !どのil行目のlatlonもilatilonと一致しなかったとき、データは存在しないため、欠損値を入れる goto 100 100 continue enddo !lon enddo !lat do mon=1,nmon write(20,rec=mon)((snow(lon,lat,mon),lon=1,nx),lat=1,ny) !rec=時間とするとgradsで読みやすい write(6,*) ((lon,lat,snow(lon,lat,mon),lon=1,nx),lat=1,ny) !結果の画面表示 enddo !mon end g95 csv-bin.f90 ./a.out 6 binaryを読むctlファイルの作成 dset ^test.bin title test options big_endian undef 1e+20 xdef 360 linear -180.00 1.00 ydef 60 linear 24.50 1.00 zdef 1 levels 1000

tdef 12 linear jan1986 1mo vars 1 q 0 99 x,y q endvars エンディアンはlittleかbigか確認する vi test.ctl grads open test.ctl d q _{で図に描ける。} 7 grads で左の図が描ける。 open test.ctl

set display color white c

set lat 50 90

color 0 100 10 -kind white->aqua->dodgerblue->black

set gxout grfill map nps d q cbar 8 Fortran について参考になるページ に掲載のリンク参照 http://seesaawiki.jp/ykamae_linux-note/d/Fortran