気象系技術メモ

気象データ解析のための技術まわりのメモ。

GrADsのコントロール（ctl）ファイルについて

GrADsで格子点データを描画する際には、データそのものを保存しておくバイナリファイルと、それがどのようなデータかを記しておくコントロールファイル（ctlファイル）というものが必要です。

自分でデータを作る場合にはこのコントロールファイルも自作する必要があるので、以下その中身についてまとめます。

なお、公式ドキュメントでctlファイルの仕様は網羅されていますので、詳しく知りたい方はこちらを参照してください。
cola.gmu.edu

たくさん設定する項目がありますが、そんなに使わないものも多いので、よく使うところをかいつまんで紹介しようと思います。



コントロールファイルの中身はこんな感じです。

DSET hoge.bin
UNDEF 9.96921e+36
OPTIONS yrev
XDEF 288 LINEAR 0 1.25
YDEF 145 LINEAR -90 1.25
ZDEF 1 LEVELS 500
TDEF 20 LINEAR 12Z03JUL2020 6hr
VARS 2
var1_name 0 99 var1_description
var2_name 9 99 var2_description
ENDVARS

それぞれの項目が何を示しているのかは以下の通り。

DSET

データが入っているバイナリファイルのパスです。

UNDEF

欠損値（undef）として処理する値を指定します。
この場合、ファイルを読み込んで9.96921e+36という値があれば、その点では値がないとして処理されます。

OPTIONS

ファイル読み込み時のオプション。
いろいろあるので、詳しくは上記公式ドキュメントを参照していただきたいのですが、よく使うものは
・yrev：データが北→南の方向に入っていることを示す（デフォルトでは南→北として解釈される）
・zrev：データが上層→下層の方向に入っていることを示す（デフォルトでは下→上として解釈される）
・template：DSETをテンプレ形式にし、それに従う連番のファイル群を一度に読み込む。詳しくはこちら
など。
特にyrev、zrevは間違えると厄介なので、データの方向は要確認です。

XDEF

X方向（経度方向）の格子間隔を示します。

XDEF a LINEAR b c

で、bからcごとに全部でa個のグリッドが存在することをしめしています。たとえば、

XDEF 288 LINEAR 0 1.25

は、経度0度から1.25度ごとに全部で288個のグリッドが刻まれていることを意味します。

このように開始点と格子間隔で指定することもできますし、グリッドを一つずつ指定することもできます。その場合LINEARではなくLEVELSを用いて、

XDEF 12 LEVELS 0 30 60 90 120 150 180 210 240 270 300 330

のようにします。

YDEF

緯度方向のグリッド数を表します。ほぼXDEFと同様ですが、ここではLINEARとLEVELS以外にガウス格子系を使うこともできるそうです。詳しくは上記公式ドキュメントへ。

ZDEF

鉛直層数を表します。XDEFと同様です。

TDEF

時間ステップ数を表します。XDEFなどと書き方は同様ですが、開始時刻は以下のように指定します（日本語での説明が難しいので、公式ドキュメントをそのまま引用します）。

hh:mmZddmmmyyyy
where:

hh = hour (two digit integer)
mm = minute (two digit integer)
dd = day (one or two digit integer)
mmm = 3-character month
yyyy = year (may be a two or four digit integer; 2 digits implies a year between 1950 and 2049)

mmはなくてもよいみたいです。

つまり、たとえば2020年7月3日12Zを初期時刻とする場合、

21Z03JUL2020

となるということです。

VARS

データファイル中に格納されている変数を記述します。まず、

VARS n

でn個の変数が保存されていることを指定し、次に改行してそれぞれの変数について記述します。

var1_name 0 99 var1_description
var2_name 9 99 var2_description
...

並び順は

varname levs units description

という形で、それぞれ
・varname：変数名（15字以内）。描画時に d varname のように指定するのに用いる。
・levs ：その変数の値が入っている層の数。地表面データであれば0とする。
・units：ファイル中の変数を指定するパラメータ。GRIB等を読み込む際に必要で、バイナリを読む際は"99"としておく。
・description：変数の説明（140字以内）。正式名称とか単位とか。

変数がすべて書き終わったら、必ず最後にENDVARSをつけます。これ以降は何も書きません。


コントロールファイルが完成したら、GrADsでは次のようにコントロールファイルを開き、描画することができます。openするのはもとのデータファイルではなくコントロールファイルです。

ga-> open hoge.ctl
ga-> d varname

GRIBデータの扱い②

前回の続き。
前回まとめたような手続きでバイナリダンプしたGRIBのデータをfortranで読み込んで解析します。

バイナリファイルの構造

こちらやこちらに詳しくまとまっています。

wgribでGRIBをバイナリダンプしたファイルには、


時間--変数--高度--緯度--経度

の順で辞書式に値が入っています。複数変数のデータがある場合、その格納順序はもとのGRIBファイルにおけるそれと同じです。

個人的には、変数ごとにファイルが独立している方が何かとやりやすいのですが（複数の解析で同じデータファイルを再利用することを考えると、ファイルがまとまっていると場合によっては使わない変数まで読み込む必要が生じるため）、その場合ファイルは「経度方向グリッド数×緯度方向グリッド数×鉛直層数×時間ステップ数」の配列であることになります。

プログラムに読み込む

したがって、GRIBから作ったバイナリファイル（1変数）をFortranなどのプログラムに配列形式で読み込むとすると、

program sample
  implicit none
  integer, parameter :: nlon = 288       !経度方向グリッド数
  integer, parameter :: nlat = 145       !緯度方向グリッド数
  integer, parameter :: nlev = 37        !鉛直層数
  integer, parameter :: ntim = 1         !時間ステップ数
  real(4) :: data(nlon, nlat, nlev, ntim)

  open(11, file=ifile, access="direct",status="old", form="unformatted", convert="little_endian", recl=4*nlon*nlat*nlev*ntim)
  read(11, rec=1) data　　　!nlon*nlat*nlev*ntimの配列dataにファイルの値を入れる
  close(11)

みたいな感じになります。

探査の方式には順番探査（access="sequential"）と直接探査（access="direct"）があり、順番探査がデータの初めと終わりを示す記号（record marker）を見てその間にあるデータを読み込む方式であるのに対し、直接探査はあらかじめ定められたデータの単位長さ（record length、変数reclで指定）ごとにデータを読み込む方式です。
つまり、直接探査の場合たとえばレコード長を4バイトとするならば、4バイトのデータを一単位として、4バイトずつ読み込んでいくことになります。
ひとつひとつの数値のデータサイズがすでに決まっているバイナリファイルでは直接探査の方が便利なので、ここではaccess="direct"としています。


上記サンプルコードではrecl=4*nlon*nlat*nlev*ntimとなっていました。ここでは読み込むファイルの数値がひとつあたり4バイトで、それが「経度方向グリッド数×緯度方向グリッド数×鉛直層数×時間ステップ数」の数だけ存在するので、データのサイズは


4バイト×経度方向グリッド数×緯度方向グリッド数×鉛直層数×時間ステップ数

ということになります。したがって、recl=4*nlon*nlat*nlev*ntimというのは「レコード長＝ファイル全体のサイズ」ということで、ファイルをまるごと一気に読み込んでしまうということです。


そのあとのread文のrec変数は「レコード番号」で、レコード長ごとに頭からデータを数えていったときのデータの番号です。今回はファイルをまるごと一つのデータ（配列）として読み込んでいるのでレコード番号は1しかありませんが、例えば「各時間の500hPa面の値だけ取り出したい」ような場合には、レコード長を


recl=4バイト×経度方向グリッド数×緯度方向グリッド数

として、

real(4) :: data(nlon, nlat, 1, ntim)  !層数=1

  open(11, file=ifile, access="direct",status="old", form="unformatted", convert="little_endian", recl=4*nlon*nlat)
  do t=1,ntim
  read(11, rec=(t-1)*nlev+22) data　　
  close(11)

とすることもできます。22というのが、鉛直層の格納順序で見た時の500hPaの順番（22番目）です（使うデータによっても異なるので要確認）。
この場合、レコード長＝レコード単位を水平グリッド数としていますから、時刻tにおける500hPa高度（22番目）のグリッドデータのレコード番号は(t-1)*nlev+22ということになります。


あと地味に厄介なのがエンディアンの指定で、リトルエンディアンとビッグエンディアンを間違えると値が発散したりしてしまうため、計算実行時には注意が必要です。逆にいえば、計算結果が明らかにおかしいときには計算ミスかエンディアンのミスが疑われます。

なお、前回紹介した方法でwgribを使ってバイナリに変換した場合、変換時のエンディアンはCPUに依存するそうです。そのため、ほかの人からもらったデータではエンディアンが違う場合があります。
用いるファイルのエンディアンがすべて共通している場合には、コンパイル時に

gfortran -fconvert=little(または big)-endian -frecord-marker=4 hoge.f90

のようにオプションをつけてまとめて指定してしまってもよいですが、ファイルごとにエンディアンが異なる場合は読み込み時のopen文で

open(11, file=ifile, access="direct",status="old", form="unformatted", convert="little_endian", recl=4*nlon*nlat*nlev*ntim)

のようにconvert="little_endian"などと個別に設定してあげる必要があります。


このようにしてデータを読み込めたら、あとは配列同士で適当に計算して欲しい値を作ってあげればOKです。


出力も読み込みと同じく、

open(12, file=ofile, form="unformatted", status="replace", access="direct", convert="little_endian", recl=4*nlon*nlat*nlev*ntim)
write(12, rec=1) output_array
close(12)

のように行います。



上記のようにして作成したデータをGrADsで描画するためには、データを書き込んだファイルとは別に「コントロールファイル」というものが必要です。
詳しくは次の記事をご覧ください。

GRIBデータの扱い①

GRIBファイルをGrADsで描画できるようにしたり、プログラムで解析したりできるようにするための操作について、忘れないうちにメモ。

手順ごとの詳細な説明というよりは、それぞれの部分で詳しいサイトに飛べるようなまとめみたいな感じかも。Mac、Unix、Linux向け。

そもそもGRIBってなによ

まあwikiを見てもらえばいいんですが、簡単に言えばWMO（世界気象機関）が導入した気象データの格納形式で、格子点上での値をバイナリで記録する規格です。現状ではGRIB1とその後継であるGRIB2の2バージョンが使われています。

よく見るところとしては、JRA-55やERA-Interimといった再解析データがGRIBファイル（.grb、.grib）で提供されています。

前提となる環境構築

こちらのページにて、気象データの解析に必要な環境構築の手順についてまとめてくださっています（Mac向け）。
qiita.com
以下、こちらで紹介されている環境（gccとか）で使うことを前提にまとめます。

GRIBファイルの変換

GRIBはそのままでは描画・解析ができない/難しい場合が多いので（対応しているソフトウェアももちろんある）、何らかの手順で別のファイル形式に変換する必要があります。

① netCDFへ

netCDFは、アメリカの研究機関（UCAR Unidata）により導入された汎用配列データ形式です（拡張子.nc）。

netCDFに変換することで、GrADsでかなり簡単に描画できるようになります。落としてきたデータを解析前にサクッと確認したいときや、とりあえず何か描いてみたいときにおすすめ。

変換にはCDOと呼ばれるソフトウェアを用います。CDOはコマンドラインで動作し、netCDF/GRIBの変換のほか格子系の変換や統計解析などかなり色々なことができるツールです。詳しくはこちらのページなどでまとめてくださっています。

www.hysk.sakura.ne.jp

インストール手順はこちら（上記ページ中のリンクから飛べます）で紹介されています。Macの方は公式ページに掲載されているMacPortsやHomebrewを使ったインストール方法も簡単かと思います。


CDOを使ったGRIBの変換は以下のコマンドで行います。

$ cdo -f nc copy [入力（GRIB）ファイル名]　[出力（netCDF）ファイル名]


そして、GrADsでの描画コマンドは以下のようになります。

ga-> sdfopen hoge.nc   #netCDFファイルをオープン
ga-> d 変数名           #描画

②バイナリファイルへ

単純なバイナリファイルにする方法。プログラムに読み込んで解析する際に扱いやすい形式です。

ここではwgribというツールを使います。wgribには無印のwgribとwgrib2があり、それぞれGRIB、GRIB2のファイルに対応しています（wgrib2はGRIBにも対応しているとのことですが、あまりうまくいかないことが多いようです）。ここではwgribを使ったGRIBファイルの操作について書きます。

wgribはNOAAのサイトから落としてくることができます。インストール手順についてはこちらのページがわかりやすいと思います。

tarファイル解凍→ソースコード生成→make という簡単な手続きです。wgribという実行ファイルが作られるはずなので、あとはそれを/usr/local/binなど適当な（PATHが通っている）ところにコピーすれば完了です。

wgribも先ほどのCDO同様コマンドラインで動作するソフトウェアです。wgribのコマンドはこちらのサイトが非常に参考になります。
www.hysk.sakura.ne.jp
例えば、あるGRIBファイルから可降水量のデータを切り出して、バイナリファイルに書き出したいとします。GrADsやプログラムでも読めるようにヘッダなしのバイナリファイルに出力する（ヘッダがあると、頭に余分な数値が読み込まれてしまう）には、以下のようなコマンドを入力します。

wgrib [入力ファイル名] | grep ":PWAT:" |  wgrib  -bin -nh -i [入力ファイル名]  -o [出力ファイル名]

ここで"PWAT"というのが可降水量を表す変数名で、

wgrib [入力ファイル名]

とすると、

wgrib anl_column125.2020070100 

Undefined parameter table (center 34-241 table 200), using NCEP-opn
1:0:d=20070100:PWAT:kpds5=54:kpds6=200:kpds7=0:TR=0:P1=0:P2=0:TimeU=1:atmos col:anl:NAve=0
2:62748:d=20070100:COVTM:kpds5=152:kpds6=200:kpds7=0:TR=0:P1=0:P2=0:TimeU=1:atmos col:anl:NAve=0
3:125496:d=20070100:CAPE:kpds5=157:kpds6=200:kpds7=0:TR=0:P1=0:P2=0:TimeU=1:atmos col:anl:NAve=0
4:188244:d=20070100:HLCY:kpds5=190:kpds6=200:kpds7=0:TR=0:P1=0:P2=0:TimeU=1:atmos col:anl:NAve=0
5:250992:d=20070100:PROB:kpds5=191:kpds6=200:kpds7=0:TR=0:P1=0:P2=0:TimeU=1:atmos col:anl:NAve=0

こんな感じでそのGRIBファイルに含まれる変数リストを表示してくれます。

どれがなんの変数かわからない場合は、

wgrib -V [入力ファイル名]

とすれば、

wgrib -V anl_column125.2020070100 

Undefined parameter table (center 34-241 table 200), using NCEP-opn
rec 1:0:date 2020070100 PWAT kpds5=54 kpds6=200 kpds7=0 levels=(0,0) grid=255 atmos col anl:
  PWAT=Precipitable water [kg/m^2]
  timerange 0 P1 0 P2 0 TimeU 1  nx 288 ny 145 GDS grid 0 num_in_ave 0 missing 0
  center 34 subcenter 241 process 201 Table 200 scan: WE:NS winds(N/S) 
  latlon: lat  90.000000 to -90.000000 by 1.250000  nxny 41760
          long 0.000000 to -1.250000 by 1.250000, (288 x 145) scan 0 mode 128 bdsgrid 1
  min/max data 0.0899731 69.4337  num bits 12  BDS_Ref 0.0899731  DecScale 0 BinScale -5

rec 2:62748:date 2020070100 COVTM kpds5=152 kpds6=200 kpds7=0 levels=(0,0) grid=255 atmos col anl:
  COVTM=Covariance between v and T [K*m/s]
・
・
・

のように詳細な説明と一緒に表示してくれます。

ヘッダをつけない場合、出力されたバイナリファイルのサイズが数値のバイト数*東西グリッド数*南北グリッド数*鉛直グリッド数になっていれば問題ありません（たとえば、数値が4バイト実数で格納されていれば「4*東西*南北*鉛直」のファイルサイズになる）。ヘッダありの場合、ここに先頭のヘッダ分のバイトが加わります。

注意点として、wgribはデフォルトでは北→南の方向にデータを出力するので、南→北の順になっているデータと合わせて解析する場合などは適宜（解析前にひっくり返すとかして）対応が必要です。
厄介なのが、GRIB2を扱うwgrib2ではデフォルトで南→北になっていて、北→南で出力するためには別途オプションで「-order we:ns」とつけなければならない点。
不安があれば解析前に単体で描画してしまうとミスに気付きやすいのですが、例えばGrADsのctlファイル（別の記事でまとめてます）にもグリッドの南北、上下を反転させて読み込むオプションがあったりするので、案外面倒です。wgrib用のスクリプトやctlファイルを流用する際には注意しましょう。



以上、まずはGRIBファイルを扱いやすい形に変換する手順についてまとめました。これらを使った解析については、長くなるので次の記事でまとめます。