[気象衛星 ひまわり 8号のデータ (千葉大による緯度経度格子データ) を読む] で、 気象衛星 ひまわり8号 のデータ (千葉大学 環境リモートセンシング研究センター (CEReS) で緯度経度格子に変換されたもの) を読んで作図することについて、説明した。
千葉大 CEReS では、過去の、ひまわり 1～7号や、 ひまわり 5号 (GMS-5) と 6号 (MTSAT-1R) の間の期間、 アメリカ合衆国から借りて運用されていた GOES 9号の画像データも アーカイブして提供している。
GOES 9号のデータには、上記の「静止気象衛星ひまわりアーカイブ・ダウンロード情報のページ」の 「静止気象衛星ひまわり5号およびGOES9号（GMS-5 SVISSR と GOES-9 Product）情報」 のところに 説明がある。そのドキュメントには、 データが置かれている FTPサイトは ftp://gms.cr.chiba-u.ac.jp/GOES9-Pacific/CEReS_products/ と書かれている。
実際に anonymous FTP で ftp://gms.cr.chiba-u.ac.jp にアクセスしてみた。 (立正大学ネットワーク内からはとくに許可をうけたマシンからでないと FTP がつかえない。 わたしは自宅からアクセスした。) 「GOES9-Pacific/CEReS_products/」の下にドキュメントどおりのデータがあるが、 その横ならびの 「GOES9-Pacific/gridded」の下に、ひまわり8号のばあいとほぼ同様に 緯度経度格子に変換されたデータ (タイムスタンプは 2013年) があるので、そちらをつかうことにした。
gridded の下には、たとえば「200501」のように名まえが 年・月の数字でできているディレクトリがあり、 その中にあるファイルは、たとえば 「200501312323.vis.tar.bz2」 「200501312323.ir.tar.bz2」のように名まえがついている。 数字の部分は、年、月、日、時、分 (UTC) である。 だいたい1時間ごとに vis, ir それぞれのファイルがある。
卒業研究で新潟県の大雪をあつかっている学生が 2005年 2月 1日 に注目しているので、 200501311423 から 200502020023 までのファイルをダウンロードした。 日本時間でいうと 1月31日 23:23 から 2月2日 09:23 である。
ここからの作業は、Linux 上でおこなった。
まず、ダウンロードしたファイルの名まえのリストをつくった。(ls -1 の出力を redirect した。)
「tar.bz2」の形式のアーカイブファイルは、tar と bunzip2 がインストールされていることを前提として、 「tar jxvf アーカイブファイル名」で展開することができる。 ダウンロードされたファイルのそれぞれについてこのコマンドを実行するシェルスクリプトを (awk をつかって) つくり、実行した。
irのほうのアーカイブからは、たとえば 「g09.200501311423.ir?.geoss」の「?」のところが 1, 2, 3 になった3つのバイナリデータファイルと 「CalibTable.txt」というテキストファイルが出てくる。 バイナリデータファイルの大きさは 9000000 バイトである。
visのほうのアーカイブからは、たとえば 「g09.200501311423.vis.geoss」のような名まえのバイナリデータファイルと 「CalibTable.txt」というテキストファイルが出てくる。 バイナリデータファイルの大きさは 144000000 バイトである。
CalibTable.txt は、バイナリデータファイルのカウント値から物理量に換算するための表にちがいない。 全部を確認したわけではないが、見たかぎりでは、どのアーカイブファイルからとりだしたものも 内容は同じである。ここにもコピーを置く。
ここでは、赤外チャネルのうちの IR1 と、可視チャネル VIS を読むことにした。
ファイル名の「分」のところは「23」であることが多いが、すべてがそうではない。 (年月日時が同じで分だけがちがうファイルはない。) つぎに GrADS で読むとき、ファイル名が規則的についているようにしたいので、 「g09.200501311423.ir1.geoss」 に 「g09.2005013114」、 「g09.200501311423.vis.geoss」 に 「visg09.2005013114」 のような別名 (symbolic link) をつける。 (「ln -s」をよびだすシェルスクリプトを (awkをつかって) つくり、実行した。) 展開は expand というディレクトリで実行し、展開結果がそこに置かれるようにした。
緯度経度格子にした ひまわり8号のデータとの類推で、 赤外は、緯度・経度それぞれ 0.02度 × 3000画素、 可視は、緯度・経度それぞれ 0.005度 × 12000画素 各画素 1バイト であると推測した。 緯度・経度範囲を適当に仮定して、GrADS の ctl ファイルを作成した。 ただし、可視画像は、GrADSにとって画素数が多すぎてエラーになり、 ファイルの前半に北半球がはいっていることがわかったので、 北半球だけをよむ ctl ファイルを作成した。 画像を表示させて、GrADSの海陸分布との比較から、緯度・経度の端の値をきめた。
↓ 赤外 (チャネル 1) 画像の ctl ファイル ir1.ctl
DSET ^../expand/g09.%y4%m2%d2%h2 UNDEF 255 options yrev template TITLE GOES9 count values XDEF 3000 LINEAR 80.02 0.04 YDEF 3000 LINEAR -59.98 0.04 ZDEF 1 linear 1 1 TDEF 35 LINEAR 14z31jan2005 1hr VARS 1 count 0 -1,40,1 GOES9 count value endvars
↓ 可視画像 (北半球) の ctl ファイル vis_nh.ctl
DSET ^../expand/visg09.%y4%m2%d2%h2 UNDEF 255 options yrev template TITLE GOES9 count values XDEF 12000 LINEAR 80.005 0.01 YDEF 6000 LINEAR 0.005 0.01 ZDEF 1 linear 1 1 TDEF 35 LINEAR 14z31jan2005 1hr VARS 1 count 0 -1,40,1 GOES9 count value endvars
ここでは日本とその付近の総観規模の天気を見たいので、 北緯 20～50度、東経 120～150度の範囲を切り出すことにした。
GrADS で、上記の ctl ファイルを前提として、つぎのようなスクリプトを実行した。 (ディスク容量やディレクトリ構成についての注意は、 必要を感じた人からの質問に応じて説明することにする。)
結果は、各画素が 4バイトの浮動小数点数のバイナリデータとして、 たとえば 「ir1_japan/ir1_japan2005013114.r4l」(2250000 バイト)、 「vis_japan/vis_japan2005013114.r4l」(36000000 バイト) のようなファイル名で保存される。
つぎのような ctl ファイルをつくると、切り出されたファイルを GrADS で読んで 作図することができる。
↓ ir1_japan.ctl
DSET ^../ir1_japan/ir1_japan%y4%m2%d2%h2.r4l UNDEF -9.99e8 options template TITLE GOES9 count values XDEF 750 LINEAR 120.02 0.04 YDEF 750 LINEAR 20.02 0.04 ZDEF 1 linear 1 1 TDEF 35 LINEAR 14z31jan2005 1hr VARS 1 count 0 0 GOES9 count value endvars
↓ vis_japan.ctl
DSET ^../vis_japan/vis_japan%y4%m2%d2%h2.r4l UNDEF -9.99e8 options template TITLE GOES9 count values XDEF 3000 LINEAR 120.005 0.01 YDEF 3000 LINEAR 20.005 0.01 ZDEF 1 linear 1 1 TDEF 35 LINEAR 14z31jan2005 1hr VARS 1 count 0 0 GOES9 count value endvars
カウント値は、赤外のばあいは輝度温度、可視のばあいは反射 [注] を仮にあらわすものであり、 本来の物理量への変換テーブルが、CalibTable.txt である。
【[注] 本来、反射率あるいは「アルベド」とは、反射放射の強さを入射放射の強さでわったものである。 ところが、衛星観測の可視センサーが観測された量をさすときは、 反射放射の強さを、太陽が真上にあるばあいの入射放射の強さ (一定値) でわったものを そう呼んでしまうことがある ([アルベド (albedo) (2014-07-02)] の記事で「衛星観測独特の用法」と呼んだもの)。 わたしはこれを「アルベド」(albedo) とか「反射率」(reflectance) と呼びたくない。 くるしまぎれに、単に「反射」 (reflection) と呼んでおく。 「規格化された反射」(normalized reflection) というべきかもしれない。】
日本付近を切り出した画像データと CalibTable.txt を読んで 物理量に変換された画像データを書きだすプログラムを Python で作成した。
変換結果のファイル名は、赤外は「tbb_japan/tbb_japan2005013114.r4l」、 可視は「refl_japan/refl_japan2005013114.r4l」のようにした。
変換結果のファイルを読む GrADS ctlファイルをつくった。
↓ tbb_japan.ctl
DSET ^../tbb_japan/tbb_japan%y4%m2%d2%h2.r4l UNDEF -1.0 options template TITLE GOES9 IR1 XDEF 750 LINEAR 120.02 0.04 YDEF 750 LINEAR 20.02 0.04 ZDEF 1 linear 1 1 TDEF 35 LINEAR 14z31jan2005 1hr VARS 1 tbb 0 0 GOES9 IR1 blackbody temperature [K] endvars
↓ refl_japan.ctl
DSET ^../refl_japan/refl_japan%y4%m2%d2%h2.r4l UNDEF -1.0 options template TITLE GOES9 count values XDEF 3000 LINEAR 120.005 0.01 YDEF 3000 LINEAR 20.005 0.01 ZDEF 1 linear 1 1 TDEF 35 LINEAR 14z31jan2005 1hr VARS 1 refl 0 0 normalized reflected visible radiation [non-dim] endvars
GrADSを開始し、つぎのように実行すれば、いちおう、色塗りの地図のアニメーションになる。
open tbb_japan.ctl set t 1 35 set gxout grfill d tbb
(止めるには、GrADSを起動した端末画面で control-C を押し、 GrADSプロンプトになったら quit 。)
しかし、これでは、 (色は、値が大きいところが赤寄り、小さいところが青寄りになっているが) 色と数値とのくわしい対応がわからないし、 観測日時は、右下に表示されるものの、字が小さいうえに、文字列の途中で切れてしまっている。
そこで、GrADS スクリプトをつかって、表示させることにする。
輝度温度については、 ひとまず、ひまわり 8号についてやったのと同様に、 およそ中央になる 250-260 K を薄い黄色、 そこから輝度温度が高いほうを赤系統、低いほうを青系統の段彩にしてみた。
ここでは、赤は地表面 (陸面・海面) または背の低い雲、 青は背の高い雲に、だいたい対応する。青が濃いほど雲頂が高い。 (赤外画像のふつうの白黒表示では明るいほど雲頂が高いが、ここではそうでないことに注意。)
1時刻ごとにとめたいばあいや、アニメーションにしたいばあいなどもあるが、 ここでは、1時刻の図ができるごとに画像ファイル (PNG) を書きだすことにしてみた。
結果の画像ファイルの名まえは「tbb_japan_png/tbb2005013114.png」のようにつけられる。
可視についても、ひとまず同様に、反射の値の およそ中央になる 0.25～0.3 を薄い黄色、 そこから反射が大きいほうを赤系統、小さいほうを青系統の段彩にしてみた。
結果の画像ファイルの名まえは「refl_japan_png/refl2005020100.png」のようにつけられる。 (夜はまっくらなので、朝の例をしめした。)
物理量に変換されたデータを読み、地図上に表示する Python プログラムをつくった。
IR1 輝度温度 2005年 1月 31日 14 UTC (23 JST) tbb2005013114_contourf_japan.png
VIS 反射 (無次元) 2005年 2月 1日 00 UTC (09 JST) refl2005020100_contourf_japan.png
