気象衛星 GOES 9号のデータ (千葉大による緯度経度格子データ) を読む


[気象衛星 ひまわり 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 である。

【2022年 9月に作業したときには、Linux上で、GrADS で緯度経度領域の切り出し、 Python で物理量 (赤外は輝度温度、可視は反射) への変換と作図をした。 そのときの作業記録は、このページの下の端のリンクのさきにおく。 しかし、(ひまわり8号のばあいとちがって GOES 9号については) GrADSで複数のファイルを続けて読むためには千葉大 CEReS でつけられた ファイル名をつけなおす必要があるなど、作業の流れがわかりにくかった。】

2022年11月に、(「アーカイブの展開とファイル一覧の作成」を別として、そのあとの) 一連の作業を Python でするように組みかえ、 Windows で実行することができた。ここではその流れにそって説明する。


アーカイブの展開とファイル一覧の作成

【この部分は、Linux で作業したので、その手順を書いてある。 追って、Windows での作業手順におきかえる予定。 Windows 10 のばあい、tar コマンドはあるが、bunzip2 コマンドを追加する必要がある。】

まず、ダウンロードしたファイルの名まえのリストをつくった。(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 を読むことにした。 つぎにのべる Python スクリプトを置くディレクトリの横ならび (同じ親ディレクトリの下) の 「expand」というディレクトリに、 展開された「*.ir1.geoss」の形の名まえのファイル、「*.vis.geoss」の形に名まえのファイル、 CalibTable.txt が置かれているようにする。

expand ディレクトリの中で、(Linux ならば) つぎのようにして、赤外、可視それぞれのファイル名一覧のファイルを作成する。

ls -1 *.ir1.geoss > ls1_ir1.txt
ls -1 *.vis.geoss > ls1_vis.txt

また、ファイル名から、日時 (年、月、日、時、分) の部分の文字列だけをぬきだした ファイルをつくっておく。

【ここであつかうファイル名はいずれも ASCII コードの範囲内の文字で構成されている。 Linux で上のようにしてつくれば、ファイル名一覧の文字コードも ASCII になる。 Windows でもファイル名一覧をつくることはできるが、まだ確実な手順を見つけていない。 Anaconda Powershell Prompt あるいは Windows 10 のコマンドプロンプト内で、 標準出力を redirect してテキストファイルをつくったら、 文字コードが「UTF-16」や「UTF-8」になって、つぎの動作にさしつかえた。 そのテキストファイルを、文字コードを判別できるエディタ (たとえば「サクラエディタ」) に読みこんで、 文字コード「Latin1」で保存しなおせば、つぎの作業につかうことができる。】

上の画像ファイルリストや日時文字列のリストは夜の可視画像もふくんでいるが、 夜の可視画像はほとんど情報がないので、 可視画像のリストは、あらかじめ昼間の画像をえらんでからつぎに進んだほうがよい。


緯度経度範囲を指定した切り出しと、物理量への変換 (Python による)

千葉大 CEReS で緯度経度格子に変換された GOES 9 データの 空間範囲は、南緯 60度 ~ 北緯 60度、東経 80度 ~ 200度 (西経 160度) である。 画素間隔は、赤外チャネルでは緯度・経度 0.04度 (約 4 km)、 可視チャネルでは緯度・経度 0.01度 (約 1 km) である。

ここでは日本とその付近の総観規模の天気を見たいので、 北緯 20~50度、東経 120~150度の範囲を切り出すことにした。 画素数は、赤外では 750×750、可視では 3000×3000 となる。

GOES 9号のデータは、各画素 1バイトの 0から 255 の間の整数値で収録されているが、 それは、赤外のばあいは輝度温度 (単位 K)、可視のばあいは反射 (無次元) [注] を仮にあらわすものであり、 本来の物理量への変換テーブルが、CalibTable.txt である。

ここでは、つぎの Python スクリプトによって、切り出しと物理量への変換を つづけておこなう。 できたファイルは、作業ディレクトリの横ならびの、 赤外は「tbb_japan」、可視は「refl_japan」のディレクトリに置かれる。

結果は、各画素が 4バイトの浮動小数点数のバイナリデータとして、 たとえば (2005年 2月 1日 03:23 UTC のばあい) つぎのようなファイルで保存される。


物理量に変換されたデータの作図 (Python による)

Python で、物理量に変換されたデータを読み、地図上に表示するプログラムをつくった。

緯度・経度の範囲は、つぎの2とおりをこころみた。

Python スクリプトはつぎのものを試作した。 いずれも、上でつくった日時のリストのファイルを読んで、 それにふくまれた日時についてつづけて処理するようにした。 1時刻の図ができるごとに画像ファイル (PNG) を書きだすことにした。 画像ファイル名とそれが置かれるディレクトリ名は、スクリプトを見てほしい。

物理量の階級分けと色は、まったく暫定的なもので、もっとくふうしたほうがよいだろう。

輝度温度については、 ひとまず、ひまわり 8号についてやったのと同様に、 およそ中央になる 250-260 K を薄い黄色、 そこから輝度温度が高いほうを赤系統、低いほうを青系統の段彩にしてみた。

ここでは、赤は地表面 (陸面・海面) または背の低い雲、 青は背の高い雲に、だいたい対応する。青が濃いほど雲頂が高い。 (赤外画像のふつうの白黒表示では明るいほど雲頂が高いが、ここではそうでないことに注意。)

可視についても、ひとまず同様に、反射の値の およそ中央になる 0.25~0.3 を薄い黄色、 そこから反射が大きいほうを赤系統、小さいほうを青系統の段彩にしてみた。


作図結果の例 (いずれも sjapan領域)


IR1 輝度温度 2005年 2月 1日 03:23 UTC (12:23 JST) tbb2005020103_contourf_sjapan.png

tbb2005020103_contourf_sjapan.png


VIS 反射 2005年 2月 1日 03:23 UTC (12:23 JST) refl2005020103_contourf_sjapan.png

refl2005020103_contourf_sjapan.png


画像間の演算をして作図する例: 輝度温度の時間変化 (sjapan 領域)

原則として1時間の変化を計算しようとしたのだが、 次の例は、中間の時刻で有効な画像がとれていないので、 約2時間の変化をしめしている。単位時間あたりにはしていない。 図の表題 (プログラムの変数 stitle の内容) もそれにあわせてかえている。

プログラムは、1組のデータを処理して1つの図をかくものとし、 読むファイルの日時をプログラム内の変数 sdatetimeb, sdatetimef で指定している。

IR1 2時刻の輝度温度の差 2005年 1月 31日 14:23 UTC から 16:11 UTC (表題と画像ファイル名は後の時刻でつけてある)

tbbdiff200501311611_contourf_sjapan.png

赤は輝度温度の増加で、雲が消えたり雲頂高度が低くなったところと推測される。 青は輝度温度の減少で、雲が出現したり雲頂高度が高くなったところと推測される。


旧版ページ (GrADS 併用) [気象衛星 GOES 9号のデータ を読む [2022年9月版]]


2022-11-18, 2022-12-17
増田 耕一 (MASUDA Kooiti)
[計算機を利用した気象分野の実習の教材]