京都の最近のレーダーアメダスデータの解析
最近のレーダーアメダスデータは、経度方向1/80度、緯度方向1/120度グリッドの制度を持っている。だいたい1km四方の精度ということになる。これを京都北部-琵琶湖にかけての領域について解析してみる。 以前四国周辺エリアについて行った解析プログラムをそのまま利用する。この時は、北緯32.0-35.0度、東経131.0-135.0度を切り出していた。今回は、北緯33.0-36.0度、東経134.5-137.5度を対象にしてみようと思う。
作業はデータサーバで行う。作業ディレクトリは、/hoge/terao/radar_amedas。データ展開用ディレクトリkinki_res3を作成し、領域定義ファイルdefs.txtを作成する。
type=lonlat x1=134.5 nx=320 y1=36.0 ny=360 nodata=-999.9
作業を行うためのshell script extract_kinki_res3.20170107.shを作成してこれを実行する。dvdのisoイメージを各年ごとに2006-2013年についてマウントし、作業を行う。その結果を平均することで、kinki_rese_productsディレクトリに様々な平均ファイルを作成する。
2006年1月2日05JSTのデータのみが欠けている。
また、2011年のデータはどうも形式が違うようで、異なるscriptを使う必要がある。
これを実際に描画する作業は、解析サーバで行う。作業ディレクトリはtopics/kinki_climate/analである。ここにkinki_res3_productsの内容をコピーしてくるとよいだろう。 まずは、全期間平均と月平均を作成。 季節を決めて日変化も作成してみたいところだ。
標高データはucsdからSRTM30_PLUSをとってくる。
http://topex.ucsd.edu/cgi-bin/get_srtm30.cgi
上記の解析エリアより少し広い領域を採用。
データの描画はGMTで行う。 降水量のグリッドデータはただの4バイトバイナリである。緯度北から南、経度西から東で入っている。 データ数が360x320個、レコード長460800バイトのファイルとなっている。 地形のデータは30秒ごと3カラムのテキストデータで、各カラムが経度(deg)、緯度(deg)、標高(m)となっている。
まずは、降水量データをgrdimageで書く。
その上に、grdcontourで地形のコンターを重ねつつ、海岸線と緯度経度を置くのが良いだろう。
grdimageとgrdcontourで許されるグリッドデータの形式はどのようなものがあるか、GMTのマニュアルを参照する。
http://gmt.soest.hawaii.edu/doc/latest/index.html
Cのnativeなバイナリもサポートしているようだから、単純な4バイト浮動小数点実数も大丈夫だと思われる。
ただ、データはgrd形式のファイルが便利そうではある。
地形データ(x,y,z)からはどうやらxyz2grdコマンドによってgrd形式のファイルをつくれそうである。
単純な浮動小数点実数データからも、xyz2grdコマンドのオプション指定により、grd形式のファイルが作れそうである。
まずは以上の見通しをもって、作業にあたった。
枠と海岸線を書く
まず最も簡単そうなことから。枠と海岸線を書いてみる。 参考にしたのはまずはこのページ。
http://www.creator.club.ne.jp/~kami/tips/gmt/japan.html
マニュアルはフルにあらゆる使い方に対応するコマンドラインオプションを説明してくれているため参照するにはちょっと重すぎ。 ふつうこの程度で使うのだぞ、という情報としてはこのページの方が分かりやすい。以下を実行してみる。
gmt pscoast -R134.5/137.5/33.0/36.0 -JM12c -Df -W2 -Ba1/f0.5/g0.2 > kinki.eps
まずはまったのは、
ペンの指定の-Wオプション。
上記ページの例示には、-W1/0/0/250という記述がなされている。
ところが、これは正しく解釈されないという警告が出る。おそらく書式が古い推奨されないのだろうと思う。
まずは-W2で動くかどうかやってみることにしたら、いちおう警告は消えた。
次にはまったのは、
pscoastコマンド(最近はgmt pscoastとラッパーに包んで使う)が海岸線データを見つけられない。
という事態。
上記-Dfオプションが、full resolutionの海岸線データを使うという意味。
ほかにも、h, i, l, cがあり、この順で粗くなっていく。
とりあえず、i ならば対応できるようだ。
エラーメッセージの中に GSHHG がないという表示が出るので、おそらく解像度の高い海岸線データがインストールされていないのだろうと思い、gshhgをyum searchしてみたところ、
[terao@gkmodel anal]$ sudo yum search GSHHG gshhg-gmt-nc4.noarch : Global Self-consistent Hierarchical High-resolution
ところが、
[terao@gkmodel anal]$ sudo yum install gshhg-gmt-nc4 パッケージ gshhg-gmt-nc4-2.3.3-1.el7.noarch はインストール済みか最新バージョンです
だから既にインストールされているが、gmtが見つけられなくなっている、ということのようだ。 環境変数を指定する必要があるのだろうと想像する。
ちなみに-Di指定で見ると、この程度。
いまいち・・・ですね。 そこで、あくまでも full resolution を目指すことに。 一体問題のファイルはどこに入っているのだろう。 パッケージの中にあるファイルのインストール位置を知るコマンドがあったはずだ・・ いろいろ検索してようやく思い出す。
[terao@gkmodel anal]$ repoquery --list gshhg-gmt-nc4 /usr/share/doc/gshhg-gmt-nc4-2.3.3/README.TXT /usr/share/gshhg-gmt-nc4 /usr/share/gshhg-gmt-nc4/binned_border_c.nc /usr/share/gshhg-gmt-nc4/binned_border_i.nc /usr/share/gshhg-gmt-nc4/binned_border_l.nc
ありゃ。full resolutionはないぞ。
[terao@gkmodel anal]$ yum search gshhg gshhg-gmt-nc4-full.noarch : GSHHG - full resolution gshhg-gmt-nc4-high.noarch : GSHHG - high resolution
あ、そういうこと。解決。
地形データのgrdファイル化
地形データはテキストファイルなので、xyz2grdサブコマンドでgrdファイル化することにする。
[terao@gkmodel anal]$ gmt xyz2grd topo/shikoku.txt -Gtopo/shikoku.grd xyz2grd: Syntax error: Must specify -R option xyz2grd: Syntax error -I option: Must specify positive increment(s)
とのことなので、-R, -Iオプションを検討する。
-Iは、gtopo30を持ってきているのだから、30秒(30 arcseconds)に設定するべきだろう。つまり-I30s/30s。
-Rは、今書きたいと思っている領域にするなら、-R134.5/137.5/33.0/36.0が良いだろう。
[terao@gkmodel anal]$ gmt xyz2grd topo/shikoku.txt -Gtopo/shikoku.grd -I30s/30s -R134.5/137.5/33.0/36.0
これで地形のgrdファイルtopo/shikoku.grdができた模様。
地形をコンターで描く
コンターを書くならgrdcontourである。
gmt pscoast -R134.5/137.5/33/36 -JM12c -Df -W2 -Ba1f0.5g0.2 -K > kinki_topo.eps gmt grdcontour topo/kinki.grd -R -JM -W1 -C200 -L200/3000 -O >> kinki_topo.eps
うまくいくかと思いきや、妙なことになってしまった。
たぶんgrdファイルが特定の格子点上のグリッドだけについて変換されたのではないか、と推察する。 これらのグリッドだけコンターがかかれ、他がゼロだとこうなるだろう。
グリッドの定義は通常gridline registrationであることと関係がありそうだ。 topo30は、pixel registrationでデータが定義されている。
gridline registrationで、かつグリッド幅を、-Rオプションで指定された領域の端がグリッドになるように調整すると、もしかすると等間隔にのみグリッド地が定義されるような状況にもなりうるだろう。
xyz2grdに-rオプションを指定して実行してみる価値があるだろう。
やってみた。
[terao@localhost anal]$ gmt xyz2grd topo/kinki.txt -Gtopo/kinki.grd -I30s/30s -R134.5/137.5/33.0/36.0 -r [terao@localhost anal]$ gmt pscoast -R134.5/137.5/33/36 -JM12c -Df -W0.5 -Ba1f0.5g0.2 -K > kinki_topo.eps [terao@localhost anal]$ gmt grdcontour topo/kinki.grd -R -JM -W0.3 -C200 -L200/3000 -O >> kinki_topo.eps
見事解決。
いよいよデータの描画
データファイルは4byte floatの塊になっている。
これをgrd形式のファイルに持ち込めばあとは簡単である。
そのためには、xyz2grdが使える。具体的には、
[terao@localhost anal]$ gmt xyz2grd rain.m08.bbx -Grain.m08.grd -I0.0125/30s -R134.5/138.5/33.0/36.0 -r -ZTLf
のような感じである。
-Iオプションはグリッド間隔。
-Rオプションはグリッド定義で、-rオプションでpixel registrationを指定。
-Zオプションは、データがどういう向きに入っているかを示す。
この場合、Tつまり、上=北から南。
さらにその各行が、Lつまり、左=西から東。
なお、作成している平均値ファイルは4レコードから構成されているため、
上記xyz2grdだと「余分なデータがある」とはじかれてしまう。
そこで、必要な部分だけを切り出すため、headコマンドを使った。
headコマンドは、-c Nオプションをつけることで、
初めからNバイトの部分を切り出してくれる。
今回の場合、320 x 360グリッドなので、320 x 360 x 4=460,800バイトの切り出し。
すなわち、
head -c 460800
を実行すればよい。
このgrd形式のファイルに基づいてイメージを描く。grdimageを利用する。
さて、その場合、カラーパレットの指定が少し面倒である。
/usr/share/gmt/cptディレクトリに、標準的な多くのカラーパレットが定義されている。
今回は、オンライマニュアルのcook bookの22節, "Color Space: The Final Frontier" を参照してなんとなくよさそうだった、wysiwyg.cptを使ってみた。
雨の表現として結構自然だと思う。
問題はこのカラーパレットを現実のデータに対応させることだ。
このためのツールgrd2cptも用意されていて、そうとう便利。
[terao@localhost anal]$ gmt grd2cpt rain.grd -Cwysiwyg.cpt -L0/10 > hoge.cpt
上記、rain.grdを参照して、0-10mm/hの範囲で、各色の幅が等頻度で現れるように調整したカラーパレットを、hoge.cptに作成してくれる。
実際にはこれをもとにして、何らかの「切りの良い」境界値を決めてやるのがよいだろう。
まずはこれでとりあえずの作業はできるようになったように思う。
明日少し整理して、結果をお示ししたいと思う。 もう寝る時間です。
PostScriptファイルをokularで表示
あと、意外とはまったのは、psファイルはどのコマンドで見たらよいのか、という点。
以前他のマシンでpdfファイルの表示ではまったことがあったのだが、コマンドの固有名詞を忘れてしまうとどうにもならない。
いろいろ検索してようやく思い出したのがokular。
ところが、このマシンにはこれはまだインストールされていなかったため、yum installする必要があった。
153個もパッケージをインストールする必要があったため少し時間を要し、ようやく終了。
無事表示ができた。
convertコマンドでトリミング
convertは周囲の余白をトリミングしてくれる。こんな感じに使う。
[terao@gkmodel anal]$ convert -rotate 90 -resize 50% -trim kinki.eps kinki.png
備忘のため。
